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down  on  computer  time,    a   search  was  made  up  over  only  the  two  variables 
v,w,    since  u,    the  heading  angle,    is  of   relatively  small   significance  in 
the  computation  of  the  Hamiltonian.      The  search  consisted  of   taking  all 
possible  combinations  of  v,w  where  v  ranged  over  the  set  of   values   3,6, 
9,12,15  and  w  ranged  over   the  values   200,400,600,800,1000  and  comparing 
the  Hamiltonian  for  each  set  with  the  one  generated   in  the  numerical   rou- 
tine.     It    is  possible  for  the  search  with  this   size  grid  to   fail   to  de- 
tect a  set  of   controls  that  yield  a   smaller  value   for  the  Hamiltonian. 
However,   the  time  factor   is   critical    and  this   grid   size  was   felt  to   fur- 
nish a  satisfactory  compromise.      A  search  relying  on  gross   computation 
can  easily  become  completely  unrealistic    in  terms  of  the  computer  time 
required. 

The  Hamiltonian 

(3.9)  H  =  "A  ."V 

has   a  constant  tfalue  on  an  extremal,    with  allowance  being  made   for  round 
off   errors,    whenever  t   does   not  occur  explicitly   in  H. 

At   this    point,    we  might  mention  that   there   are   two   types  of   WNfEUSt* 
tions  of  the  control  variables   in  the  classical    literature,    called  weak 
variations  and  strong  variations,    [k  ~j    Weak  variations  are  variations   in 
which  the     |  </uL      are    "small"   for   each  time   steaHpstrong  variations   are 
variations   in  which  C      |</uL|    dt   is   "small".      That    is,    in  weak  variations 
only  values  of  control   near  those  used  are  compared  but   if   strong  varia- 
tions  are  considered,    then  the  new  control    function  may  not   be   "near"  the 
one  used. 

All  methods  of  determining  the  routes   are  methods  of  variations  which 
deform  a  given  path.      A  path  which  furnishes   a  relative  minimum,    if  we 
allow  only  wel%  variations ,   may  not   furnish  such  a  minimum  if   strong  var- 
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iations   are  allowed,    as  will  be  seen  in  section  6. 
5.      The  Numerical  Routine. 

The  routine  for  determining  the  route  is   given  in  this   section. 

Heuristic  Discussion: 

Let   us  guess  a  set  of  values   for  the  parameters  h   ,   h-.      We  will 
then  use  this  set  of  values  to  determine  the  control  variables   for  each 
time  t   by  the  minimum  principle  to  determine  a  route.      The  terminal   point 
thus  generated  will,    in  general,   differ  from  the  desired  one.      By  chang- 
ing the  values  of  h, ,   h     appropriately,   this  terminal  point  will  be 
forced  toward  the  desired  point  x™,   y   . 

Mathematical   Derivation; 

First  we  see  from   (4.3)    that 

(5.1)  SA    =    A1dh]>   +    T^dty 

Note  that     Va  1,    and  this   in  turn  implies  <fV  a  o.      These  facts  will  be 
used   in  the  following  equations. 

Next,    from  the  first  of  the  Euler  equations    (4.4)    by  taking  the 
total   differential  we   find  that 

-  v  sin  u/A  +  v  cos  u//'  +    (-Jv  cos   u  ~/v  sin  u  +  f     )  /u 

(5.2)  UU 

+(-A   sin  u  -  /"  cos   u  *   f     )^v  +  f     <fw  =  0, 

uv  uw 

if  we  assume  that  we  can  change  A t^  ,    u,   v,   w,   at  fixed  x,   y,    z,    t.      But 


since 


H       =   ->  v  cos   u  -/*  v  sin  u  +  f 

uu  uu 


H       =  -  A    sin  u  - y  cos  u  +  f 
uv  uv 


H       =   f 
uw         uw 


(5.2)   becomes 
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Introduction 

The  purpose  of  this  paper  was  to  develop  a  numerical  routine  for 
solving  the  submarine  routing  problem;  the  problem  is  for  the  submarine 
to  choose  a  course  that  minimizes  the  probability  that  it  will  be  de- 
tected.  Several  difficulties  arise  in  the  numerical  solution  and  sub- 
routines were  included  to  take  care  of  these. 

In  type,  the  problem  is  a  problem  of  Bolza  in  calculus  of  varia- 
tions.  It  is  complicated  by  the  fact  that  there  may  be  several  routes 
each  of  which  appear  to  be  the  solution.   They  all  begin  and  end  at  the 
desired  points  and  they  all  satisfy  the  Euler  equations.   It  is  only 
after  routines  are  incorporated  to  check  other  conditions,  the  condi- 
tions of  Legendre  and  Weierstrass,  that  it  becomes  clear  whether  a 
particular  solution  is  the  desired  solution.   When  these  conditions  are 
not  satisfied  this  fact  must  be  determined  and  a  routine  made  up  to 
determine  the  controls  to  satisfy  it. 

In  general,  the  route  is  generated  as  an  extremal,  a  solution  to 
the  Euler  equations,  though  the  basic  principle  is  that  the  control  must 
minimize  the  Hamiltonian.   The  problem  is  complicated  by  the  fact  that 
the  control  which  furnishes  the  minimum  may  lie  on  the  boundary  to  the 
region,  as  when  the  submarine  is  at  maximum  depth.   Several  subroutines 
must  be  adjoined  to  treat  the  case  when  the  control  lies  on  one  of  the 
bounding  faces  or  edges. 

The  existence  of  a  corner  introduces  further  complications;  at  a 
corner  the  control  is  a  discontinuous  function  of  time.   It  is  necessary 
to  make  up  a  search  routine  for  other  values  of  the  control  which  may 
decrease  the  Hamiltonian,  and  it  is  necessary  to  compromise  between  the 
demands  of  computing  time  and  accuracy  in  this  routine. 


I  am  indebted  to  my  advisor  Professor  F.  D.  Faulkner  for  the  en- 
couragement and  guidance  given  me.   This  is  a  continuation  of  the  study 
^l^  which  he  began.   The  principal  purpose  of  this  study  was  to  analyze 
the  numerous  difficulties  envisioned  and  to  make  up  numerical  routines 
to  resolve  these  and  others  that  developed,  so  that  a  solution  could  be 
generated  automatically  on  the  computer.   I  also  express  my  thanks  to 
Professor  W.E.  Bleick,  whose  guidelines  were  followed  in  programming 
the  numerical  routine  associated  with  this  problem,  and  to  Mrs.  Sally 
B.  Kline  for  help  when  programming  difficulties  were  encountered. 


1.  Statement  of  the  Problem. 

The  basic  problem  studied  here  is  the  following.   A  submarine 
located  at  point  xn,y  *-s  to  make  a  voyage  to  point  x  ,y  in  a  specified 
time  T.  Throughout  this  voyage  the  submarine  is  subjected  to  enemy  de- 
tection devices  whose  capabilities  are  assumed  known  statistically.   If 
the  submarine  has  previously  gone  undetected,  let  the  probability  of 
detection  in  a  time  interval  At   be  approximately  f (x,y,u,v,w,t)  At,  and 
hence  the  probability  of  being  detected  along  the  route  satisfies  the 
equation 

(1.1)  dp  =  (1  -  p)  f (x,y,u,v,w,t)  dt. 

The  function  f   is   the  best  estimate  of  the  enemy's  detection  capabili- 
ties  based  on  information  we  have  about  his   listening  devices,    tests  we 
have  run  on  our   submarines   using  comparable  devices,    the  distance   in- 
volved  in  the  trip,    and  other   information  available  to   us. 
Equation   (1.1)    can  be  simplified  by  letting 

(1.2)  z   =  -   ln(l  -  p) 
which  gives 

(1.3)  z  =   f. 

Since  we  are  primarily  interested  in  long  routes,  and  the  time  to 
change  depth,  speed,  and  heading  angle  is  assumed  small  compared  to 
total  time,  it  will  be  ignored. 

2.  Equations. 

Because  routes  of  approximately  2500  miles  are  of  primary  interest, 
the  flat  earth  assumption  will  be  used.   The  equations  governing  the 
system  may  then  be  written  as 

x  =  v  cos  u 

(2.1)  9  =  v  sin  u 

z  =  f (x,y,u,v,w,t). 


limited  range.   The  constant  c„  reflects  the  possibility  of  detection  by 
some  passing  ship,  for  example. 

The  thermocline  is  defined  as  the  depth  at  which  the  temperature 
gradient  (rate  of  decrease  of  temperature  with  increasing  depth)  is  a 
maximum.   In  equations  (2.3),  wQ  is  the  depth  of  the  thermocline,  s   is 
the  distance  from  the  terminal  point  (x  ,y  )  of  the  route  to  the  main 
concentration  of  enemy  listening  devices  which  is  represented  by  the 
point  (x2>yo)»  and  ©  is  the  angle  that  the  perpendicular  to  the  enemy 
shore  line  makes  with  the  x-axis.  The  constants  a..  ,b^  ,c,  ,d,#and  wn  are 
chosen  so  as  to  simulate,  by  the  function  f ,  conditions  as  they  exist 
in  any  given  situation. 

The  two  functions  f,  and  f  are  intended  to  be  typical  of  the 
functions  which  do  describe  the  enemy's  defenses.  The  function  f 1 , 
representing  the  passive  enemy  defense,  tends  to  decrease  as  the  dis- 
tance from  the  enemy  shoreline  increases.   This  mathematical  model  of 
the  enemy's  passive  defenses  also  tends  to  be  insensitive  in  the  region 
of  the  thermocline.   The  function  f„,  which  represents  the  enemy  active 
defense,  was  constructed  so  as  to  emphasize  the  enemy's  surface  search. 
For  this  reason,  this  function  decreases  as  the  depth  increases. 

If  submarine  routing  as  described  in  this  paper  were  to  be  made  a 
part  of  naval  operations,  generating  the  function  f  would  be  a  problem 
for  intelligence  and  engineers.   Such  things  as  the  sea  state  and  the 
nature  of  the  ocean  floor  in  the  region  where  the  submarine  is  operating 
would  then  have  to  be  included  in  the  function  f.   These  functions  would 
also  vary  with  time,  reflecting  changing  sea  state,  etc. 

The  problem  now  becomes  that  of  determining  the  control  variables 
as  functions  of  time  to  effect  the  desired  optimization.   That  is,  we 
want  to  choose  the  heading,  the  speed,  and  the  depth  so  as  to  go  from 


the  initial   point   (x_,y0)    to  the  terminal   point    (x   ,y  )    with  the  minimum 
probability  of  detection. 

By   (1.1)   the  probability  of  being  detected   along  the  route  is 

(2.5)  p(T)    =   1   -   exp  [-   z(T) 

and  hence   this   is  the  quantity  we  wish  to  minimize. 
3.      Adjoint   Equations. 

Let   us   consider   any  route  and   a  neighboring  route.      The  neighbor- 
ing route  will  be  generated  by  replacing  u,v,w  on  the  original   route 
by  u  +  <fu,   v  +  <Tv,   w  +  </w  respectively.     This   generates   first-order 
changes   in  x,y,z  satisfying  the  differential   equations 

ox  =   -  v  sin  u  cTu  +  cos   u  cfv 

(3.1)  cTy  =       v  cos   u  o u  +   sin  u  <fv 

cfz  =   f     <Ai  +   f     o v   +   f     <A>7  +   f     cTx  +  f     <Ty. 
u  v  w  x  y 

3  f 
The  notation  f  * ,  etc.  as  used  above  will  be  employed  throughout 

the  remainder  of  the  paper. 

Let  us  now  introduce  three  Lagrange  multipliers  ~\,y,S   which  are 

some  unspecified  function  of  time  t.   Multiply  equations  (3.1)  by  A,y  , 

V  in  that  order,  add  the  three  equations,  and  integrate  the  result  over 

the  interval  (0,T).   These  operations  yield 

np  — 

(3.2)  f  >(<fx  +  v  sin  u  <Tu  -  cos   u  Sv  +/'(<fy  - 

v  cos   u  ^u  -    sin  u  ^v)   +  V(^z  -   f  <Px  -   f  cfy 

x  y 


-   f     cfu  -   f  <fv  -   f  dw)        dt. 

U  V  w  J 


Separating  the  terms  containing  the  variations  of  the  state  variables 
from  those  containing  the  variations  of  the  control  variables,  we  get 
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(3.3) 


(       U  <fx  +  j-Sy  +  V(<fz  -   fx<Tx  -   f  <fy)~]  dt 

~  \      ^  |~^v  s*-n  u  +>•  v  cos   u  +  Vf    1  «Ai  +    J  A  cos  u 
Jo 


+  >fsin  u  +  Vf  \<Tv   +  Vf  </w)  dt. 
vj       w 

Integrating  by  parts  in  (3.3)  just  the  terms  on  the  left  involving 

derivatives  of  state  variables,  we  get 


(3.4) 


[*<fx  +j-<Ty  +VcTz]       -     (      f</x(> 


A  +  Vf  ) 
x 


+  /y(^  +Vf   )    +  (fz-J 

y  J 


dt 
y 


which  combined  with  the  right-hand  equations  of    (3.3)   yields 

T  XT 

I   *<Tx  +^<fy  +V<TZJ        -     \        <Tx(a  +Vf   )    + 


0  0 


•T 


•  dt    =   \      [(->  v 


(3.5)  dy(J*  +Vf  )    +  <rzV  dt   =   \      |  (-*  v  sin  u  +-Kv  cos  u 

'o 

+Vf   )  <Ai  +   (3  cos  u  +rsin  u  +Vf  )  <Tv  +  Vf  <fw     |dt. 
u  v  w        J 

Now,    choose  A  ,>>,"/ in  such  a  manner  that  they  satisfy  the  differential 

equations 

*  =  -Vf 
x 

(3.6)  ^    s    -Vf 

v-=  0. 

Equations  (3.6)  are  called  the  adjoint  equations  of  the  variation- 
al equations  (3.1).   With  this  choice  of  ^  ,^  ,V  equations  (3.1)  reduce 
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to 


(3.7) 


a  <Tx  +  y  <ry  +  v  fz 


-V£u  ]  fn  dt  ♦  \    [. 

-'o 


T        rT 

]  •  i 


>  V  sin  u  -/'v  cos   u 


^  cos  u  +  y  sin  u  +   Vf 


v]/v 


dt   + 


t 


Vf     <TW  dt. 
w 


Note  that  this  formula  gives  us  a  relation  between  the  terminal  values 
of  x,y,z  and  an  integral  made  up  of  the  variations  of  the  control  varia- 
bles ofu,  </v,  </w.   Note  also  the  important  fact  that  the  values  of 
cfx,     Jy,  <Tz   interior  to  the  interval  (0,T)  are  not  needed. 

For  convenience  in  the  sections  to  follow,  the  vector  notation 


(3.8) 


rv  cos  u 
\/  =  I  v  sin  u 


will  be  used.   The  scalar  product  of  these  two  vectors 


(3.9)     H  =  A  *  V 
is  called  the  Hamiltonian. 

Also  for  future  use,  let  us  define  the  following  three  particular 
solutions  to  the  adjoint  equations 


(3.10) 


=   1 
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4.   Conditions  for  a  Minimum. 

In  this  section,  the  necessary  condition  for  a  minimum  will  be 
given. 

If  a  path  is  to  provide  the  desired  minimum,  it  must  first  be  admis- 
sible. 

Admissibility.   One  requirement  for  a  route  to  be  admissible  is  that 
it  begin  and  end  at  the  desired  points,  i.e.,  x(0)  =  0,  y(0)  =  0  and 
x(T)  =  x  ,  y(T)  =  y  for  some  solution  to  the  differential  equations 

X  =  V  cos  u 
(2.1)     y  =  v  sin  u 
z  =  f. 
A  further  condition  for  admissibility  is  that  the  depth  and  speed  satisfy 
the  inequalities 

m-L)  max 

0  *  v  ^  v 

max 

where  w    and  v    depend  upon  the  specific  class  of  submarine  under 
max      max 

consideration. 

A  set  of  control  variables  which  are  piecewise  continuous  and  sat- 
isfy (4.1)  are  called  allowable.  Note  that  allowability  is  a  local  con- 
straint, in  terms  of  time,  on  the  set  of  control  variables.   A  path  which 
satisfies  the  above  constraints,  (2.1),  (4.1),  is  admissible. 
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Our  problem  is  to  find  the  one  route,  among  all  admissible  routes, 
such  that  z(T)  is  a  minimum.   Minimizing  z(T)  in  turn  minimizes  p(T)  and 
this  is  our  objective. 

For  a  path  to  furnish  a  minimum,  it  must  be  admissible  and  also  sat- 
isfy the  following  conditions: 

1.  Euler  Equations. 

2.  Legendre  Condition. 

3.  Weierstrass  Condition. 

4.  Envelope  Condition. 

The  envelope  condition  was  not  investigated  and  will  not  be  treated  in 
this  paper.   The  order  of  the  conditions  as  given  above  is  used  since 
this  is  the  usual  order  in  which  they  will  be  checked  in  a  problem. 

A  point  that  should  be  kept  in  mind  throughout  is  that  on  a  path 
which  furnishes  the  desired  minimum,  the  Hamiltonian  must  be  minimized 
at  each  value  of  t.   This  is  the  Weierstrass-Pontryagin  maximum  princi- 
ple with  a  change  of  sign.   To  satisfy  this  one  criterion,  the  Euler 
equations,  the  Legendre  condition,  and  the  Weierstrass  condition  must  all 
be  satisfied. 

The  Euler  equations,  a  condition  on  the  first  derivatives,  require 
that  the  Hamiltonian  have  a  stationary  value. 

In  the  Legendre  condition  the  second  derivatives  are  checked  for  a 
minimum.   2  !  The  Euler  equations  and  the  Legendre  conditions  are  both 
local  conditions  on  the  control  variables  u,v,w. 

Finally,  in  the  Weierstrass  condition  we  compare  the  Hamiltonian 
for  the  values  of  u,v,w  used  with  the  Hamiltonian  for  all  allowable  values 
of  u,v,w,  for  all  values  of  t. 

Euler  Equations.   For  a  minimum,  the  control  variables  must  be  chosen 
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so  they  minimize 

(3.9)  H  =     A  •  "V 

as   compared  with  all   allowable  controls,    for  all   t,   0-^.t^T,    for  some 
solution  A    to  the  adjoint  system  of  equations.      Let  us  consider  our   solu- 
tion A    to   the  adjoint   equations    in  the  form 


(4.2)  A    =  hx     A  x   +  h2   X2   +  h3  A 


3 

where  h^,   h2,   h.,  are  arbitrary  constants   and    A,,    A0»    A,  are  defined 
inequations    (3.10). 

Wa  are  now  faced  with  the  problem  of  having  three  constants   in  our 
solution   A   to  the  adjoint  equations.      But  we  may  choose  one  relation 
among  the  constants  h   ,   h    ,   h   .      We  chose  h     to  be  1,    and  then 


(4.3)  A  =  hx  A      +  h2  A2  +   A3. 


This  then  leaves  us  with  the  problem  of  finding  the  other  two  constants 
so   that  x(T) ,   y(T)    will   assume  the  desired  terminal  values   x    ,   y   . 

If  the  values  of  the  control  variables   lie   inside  the  domain  of 
allowed  values,    then  the  Euler  equations 


=  H     =  -X  v  sin  u  +rv  cos   u  +Vf     =  0 


du  u  u 

(4.4)  *H  =  H     =  *cos  u  +  r  sin  u  +  Vf      =0 

Jv  v  v 

-i-5  =  H     =Vf     =  0 


3w  w  w 

are  the   first   necessary  conditions   for  the  desired  minimum.      The  Euler 
equations  when  combined  with  the  adjoint   equations   (3.6)    are  referred 
to   as  the  Euler-Lagrange  equations.      Solutions  to   the  equations  of  motion 
which  also   satisfy  the  Euler-Lagrange  equations  are  called  extremals. 

The  Euler  equations, however,  do   not   insure  that  we  effect   the  desired 
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minimum  for  H;  they  are  only  necessary  conditions.  It  is  then  necessary 
to  investigate  additional  conditions,  the  next  being  the  Legendre  condi- 
tion. 

Legendre  Condition.      This    is   an   investigation  of  the  second-degree 
terms  of  the  Taylor  expansion  of  the  Hamiltonian.      If  we  can  expand  H  in 
a  series   in  du,   dv,   dw,   valid   in   some  neighborhood  of   (0,0,0),    the  first 

necessary  condition  for  H  to  be  a  minimum   is   that  H     =H     =H     =0,    as 

u         v         w 

stated   in  equations    (4.4).      The  second  necessary  condition  is   that  the 
quadratic   form 


%       H       H      \  /  du 
uu  *  guv      uw  \  / 


(4.5) 


(du  dv  dw)        H       H       H 

vu     w     vw 


dv 


kH       H       H       /  \  dw 

\    WU       WV       WW  /     \ 

be  positive  semi-definite  at  least;   we  hope   it  will  be  positive  definite. 
This  condition,   which   is  known  as   the  Legendre  condition,  [3  J  may  be  ex- 
pressed 


H      >   0 
uu 


H       H 
uu     uv 

H       H 
vu     w 


>    0 


(4.6) 


11  TJ  TJ 

UU      uv      uw 


H       H       H 
vu     w     vw 


H       H       H 

WU       WV       WW 


>    0 


for  the  case  in  which  the  control   variables   u,v,w  all    lie  interior   to 
the  domain  of  allowed  values. 

If  one  or  more  of  the  control  variables   are  on  the  boundary,   the 
conditions  are  altered  accordingly.      Consider  the  case   in  which  w  -  w 


max 


for  some  part  of  the  route.   The  conditions  then  become 


16 


H  <  0 
w 


(4.7) 


H 

uu 

H 
uv 

H 
vu 

H 

w 

>  0 


H     < 

w 

0 

H      C 

0 

V 

H      > 
uu 

o, 

whenever  w  =  */   .   If  w  =  */  „  and  H  =0,  then  we  must  check  to  ensure 

max  max      w 

that  H   >  0.  Of  course  if  H  >  0  the  minimum  is  not  on  the  boundary  at 

that  point. 

If  v  =  v   ,  conditions  equivalent  to  those  above  will  be  used, 
max 

When  both  w  =  w    and  v  =  v   ,  the  conditions  read 
max         max 


(4.8) 


and  if  H  =  0  on  the  boundary  then  we  must  check  H  >  0,  and  similarly 

W  WW  J 

for  H  and  H  . 
v      w 

If  the  Euler  equations  are  satisfied  and  if  equations  (4.6)  are  sat- 
isfied at  some  interior  point  of  u,v,w  space,  then  these  values  yield  a 
local  minimum;  H  is  smaller  at  that  point  than  at  any  other  point  u,v,w 
in  some  neighborhood.   However,  the  point  may  not  furnish  the  minimum 
value  for  H;  it  is  necessary  in  theory  to  compare  the  value  with  the  value 
for  all  allowable  values  of  u,v,w.   This  condition,  that  the  value  chosen 
minimize  H,  is  called  the  Weierstrass  condition. 

A  comparison  must  be  made  between  the  Hamiltonian  for  the  u,v,w  gen- 
erated, H(u,v,w) ,  and  for  all  other  allowable  combinations  of  control 
variables,  say,  u  ,  v  ,  w  .   If  H(u,v,w)  >  H(u  ,  v  ,  w  ),  then  the  ^  .- 
control  variables  u,v,w  must  be  replaced  by  the  new  set  u  ,  v  ,  w  , 
which  yield  the  smaller  value  for  the  Hamiltonian. 

Unfortunately  there  is  no  easy  way  to  check  this  condition.   To  cut 
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down  on  computer  time,    a   search  was  made  up  over  only  the  two  variables 
v,w,    since  u,    the  heading  angle,    is   of   relatively   small   significance  in 
the  computation  of  the  Hamiltonian.      The  search  consisted  of   taking  all 
possible  combinations  of  v,w  where  v  ranged  over  the  set  of   values   3,6, 
9,12,15  and  w  ranged  over   the  values    200,400,600,800,1000  and  comparing 
the  Hamiltonian  for  each  set  with  the  one  generated   in  the  numerical   rou- 
tine.     It    is  possible  for  the  search  with  this   size  grid  to   fail   to  de- 
tect  a  set  of   controls   that  yield  a   smaller  value   for  the  Hamiltonian. 
However,   the  time  factor   is   critical    and  this   grid   size  was   felt  to   fur- 
nish a  satisfactory  compromise.      A  search  relying  on  gross   computation 
can  easily  become  completely  unrealistic    in  terms  of  the  computer  time 
required. 

The  Hamiltonian 

(3.9)  H  =  "A  .  "V 

has   a  constant  •tfalue  on  an  extremal,    with  allowance  being  made   for  round 
off   errors,    whenever  t   does   not  occur  explicitly   in  H. 

At   this   point,    we  might   mention   that   there   are   two   types  of   W&*&'~ 

tions  of  the  control  variables    in  the  classical    literature,    called  weak 

variations  and  strong  variations.    (~4  ]    Weak  variations  are  variations   in 

which   the      |  cTuL      are    "small"   for   each   time   stejpstrong   variations    are 

variations    in  which  \      \Ju       dt   is    "small".      That    is,    in  weak  variations 

->0 

only  values  of  control  near  those  used  are  compared  but  if  strong  varia- 
tions are  considered,  then  the  new  control  function  may  not  be  "near"  the 
one  used. 

All  methods  of  determining  the  routes  are  methods  of  variations  which 
deform  a  given  path.   A  path  which  furnishes  a  relative  minimum,  if  we 
allow  only  wel%  variations ,  may  not  furnish  such  a  minimum  if  strong  var- 
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iations  are  allowed,  as  will  be  seen  in  section  6. 
5.   The  Numerical  Routine. 

The  routine  for  determining  the  route  is  given  in  this  section. 

Heuristic  Discussion; 

Let  us  guess  a  set  of  values  for  the  parameters  h  ,  h„.   We  will 
then  use  this  set  of  values  to  determine  the  control  variables  for  each 
time  t  by  the  minimum  principle  to  determine  a  route.   The  terminal  point 
thus  generated  will,  in  general,  differ  from  the  desired  one.   By  chang- 
ing the  values  of  h, ,  h  appropriately,  this  terminal  point  will  be 
forced  toward  the  desired  point  xT,  y  . 

Mathematical  Derivation: 


(5.1)         SA    =    A1dhi  +  7\r 


First  we  see  from   (4.3)    that 

,dh2. 

Note  that  V  s  l,  and  this  in  turn  implies  <fV  =  0.   These  facts  will  be 

used  in  the  following  equations. 

Next,  from  the  first  of  the  Euler  equations  (4.4)  by  taking  the 

total  differential  we  find  that 

-  v  sin  u/A  +  v  cos  uff   +  (-Av  cos  u  -/v  sin  u  +  f  )  <fu 

(5.2)  UU 

+  (-A  sin  u  -  y   cos  u  f   f  )^v  +  f  <fw  =  0, 

uv       uw 

if  we  assume  that  we  can  change  )i  ,**  ,  u,  v,  w,  at  fixed  x,  y,  z,  t.   But 

since 

H   =  -)  v  cos  u  -/'v  sin  u  +  f 
uu  uu 

H   =  -  X    sin  u  -y   cos  u  +  f 
uv  uv 


H   =  f 
uw    uw 


(5.2)  becomes 
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(5.3)  -v  sin  u  <TA    +  v  cos   u  c/V  +  H     cfu  +  H     efv  +  H     <fw  =   0. 

uu  uv  uw 


Similarly 
(5.4) 

and  since 


cos   u  <A    +   sin  u  //  +    (-?»  sin  u  +/'cos   u  +  f      )  /u 

vu 

+   f       cf v  +   f      cfw  =  0 
w  vw 


H       =  -^Isin  u  +^cos  u  +  f 
vu  vu 

H       =   f 
w  vv 


H        ■   f 

vw         vw 

(5.4)    becomes 

(5.5)  cos   u/A  +   sin  u  />  +  H     </u  +  H     cf v  +  H     cfw  =  0. 

vu  w  vw 

The  third  Euler  equation  gives   us 

(5.6)  f     </u  +  f     cfv   +   f     </w  =  0. 

WU  WV  WW 

But   since 

f        =  H 

WU  WU 

f        =  H 

WV  WV 

f        =  H 

WW  WW 

(5.6)    becomes 

(5.7)  H     <Aj  +  H     cfv   +   H     <fw  =   0. 

wu  wv  WW 

Now  consider  the  equations    (5.3),    (5.5),    (5.7)    as   three  equations   in  the 
three  unknowns     cfu,     ov,     cfw. 

H     cfu  +  H     cTv   +  H    <^w  =  v  sin  u  <f3  -  v  cos   u  </> 
uu  uv  uw 

(5.8)  H     cfu  +  H     <fv  +  H     cfw  =   -cos   u  <f*  -   sin  u  ff 

vu  w  vw 

H     cfu  +  H     <fv   +  H     cTw  =  0. 

WU  WV  WW 
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If  the  determinant  of  coefficients  of  <Tu,   <fv,  <fvr   in  (5.8)  does  not  vanish, 
we  can  solve  for  <Ai,  <fv,  /w  by  using  Cramer's  Rule  as  follows: 
Denote  the  determinant  of  coefficients  by  D,  i.e., 


(5.9)     D= 


Huu  Huv  Huw 


H    H    H 
vu   w   vw 


H    H    H 

WU    WV    WW 


Then 


A  = 


(v  sin  u  A  -  v  cos  u  <T^)   H  H 

uv   uw 

(-  cos  u  A  -  sin  u  t/>)    H  H 

w   vw 


H    H 

WV    WW 


(5.10) 


<Tv  = 


H    (v  sin  u  <f  A  -  v  cos  u  <f/0  H 
uu  uw 


H    (-  cos  u  S*     -  sin  u  <$V  )   H 
vu  vw 


H 


wu 


H 


WW 


D 


Notice  that  (fw  could  be  found  in  the  same  manner,  but  is  not  needed 

since  it  is  not  used  in  the  numerical  routine.   Next,  since  from  (5.1) 

(f)i   =  dh  and  </>  =  dh 
1  2 

A  =   -(v  sin  u  dh.    -  v  cos  u  dtO    (H       H       -  H       H     ) 
1  2  w     ww  wv     vw 


(5.11) 


-    (cos   u  dhn    +   sin  u  dh   )    (H       H       -  H       H„  ) 
1  2  uv     ww         wv     uvr 


D 
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/v  =  (v  sin  u  dh1  -  v  cos  u  dh,,)  (Hvu  H^  -  H^  H^ 

D 

+  (cos  u  dh  +  sin  u  dh„)  (H   H   -  H   H   ) 
1  2    uu  ww    wu  uw 

D 
To  simplify  the  form  of  the  equation  let  us  set 

S11  =  T(H   H   -  H   H   )  (-  v  sin  u)  +  (H   H 

L  w  ww    wv  vw  uv  ww 

H   H  )  (-  cos  u)l  /  D 
wv  uw  J 

12 
S   =  |~(H   H   -  H   H  )  (v  cos  u)  +  (H   H 

L  w  ww    wv  vw  uv  ww 

(5.12)  H   H  )  (-  sin  u)l  /  D 

wv   uw  J 

21   r 

S   =   (H   H   -  H   H  )  (v  sin  u)  +  (H,  ,  H 

L   VU   WW     wu   VW  uu   WW 


H   H  )  (cos  u)l  /  D 
wu  uw         J 


S    =  |"(H   H   -  H   H   )  (-  v  cos  u)  +  (H   H 

L   VU   WW     wu   VW  UU   WW 

H   H  )  (sin  u)1  /  D; 
wu  uw         J 

then  the  above  equations  may  be  rewritten 

«fu  =  S11  dh  +  S12  dh 

(5.13) 

/.     21        22 
<fv  =   S   dh,  +  S   dh 
1         2 

We  now  have  the  machinery  to  derive  </x(T)  ,  </y(T).   In  deriving 

</x(T)  ,  we  use  equations  (3.5)  and  a  particular  choice  for  A   >  namely 

the  A.  =10  J  as  defined  in  equations  (3.10).   This  choice  for  A  in 


v6) 


equation   (3.7)   yields 


T        -T  ^T 

\(fx       =    \      cos   u    ePv  dt   -  I       v  sin  u    <fu  dt 


0         ~0  J0 
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which  becomes,    after  noting  that    (Tx(0)    =  0  and  substituting    /u,     <TV 
from  equations    (5.13), 

01  00  11 

(5.14)  cfx(T)   =  \     [cos   u   (S       dhx   +  S       dh  )    -  v  sin  u   (S       dh 

Jo  l 

+   S12  dh  )]      dt. 

-     (°} 

Similarly,   using  /\      =  I   1  J  ,   we  find 


2 


(5.15)  <fy(T)    =  I      [sin  u    (S21  dh     + 


22  11 

S        dh  )    +  v  cos  u   (S       dh 


0 

12  -. 

+   S        dh2)Jdt. 

Equations   (5.14)    and    (5.15)    can  be  rewritten  as 


/x(T)    =  dh-^  I       J"cos   uS        -vsinuS      "I 


dt  + 


(5.16)  J0 


C     r  22  12  1 

dh„  \    cos  u  S   -  v  sin  u  S    dt 


T 

.  f  r       21  111 

ffy(T)  =  dh  I  sin  u  S   +  v  cos  u  S   J  dt 

[  r       22  12-1 

\  sin  u  S   +  v  cos  u  S   Idt. 


If  we  set 


dh2 


(  21  11 

A11=l      (cosuS  -vsinuS      )dt 

T 

22  12 

A-,  0   =  \       (cos   u  S  -  v  sin  u  S     )    dt 


l12   =  I 

J0 
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21  11 

A        =   \      (sin   u  S    "    +  v  cos   u  S      )    dt 
21 

'0 


■i 


-T 

22  12 

A„2   =  \       (sin  u   S        +  v  cos   u  S     )    dt , 

0 


equations    (5.16)    become 

cTx(T)    =  Au   dh      +  A12   dh2 
(5.17) 

<Ty(T)   =  A21  dhx   +  A22  dh2. 

Equations    (5.17)    give  us  the  mechanism   for   a  Newton   -   Raphson   iteration 

scheme   for  correcting  h, ,    h2  ,    if  we  make  the  following   substitutions 

cfx(T)    =  xl,  -   x(T) 
(5.18) 

cfy(T)    =  yT  -  y(T)  . 

Note  that   in  equations    (5.18)     <fx(T)  ,     <fy(T)    were  equated  to   the  desired 
terminal  values  minus   those  which  were   attained.      Using   equations    (5.18), 
equations    (5.17)    become 

xT  -  x(T)    =  A,,    dh,    +  A,2  dh2 


(5.19) 

v_    -    v(r>     =    A  dh       +    /. 

22  2 


yT  -  y(T)    =  A21   dh]L    +  i 


from  which  we   are  able  to   generate  corrections   for   the  constants   h.  ,   h~  • 
Note  above  that    in  equations    (5.19)    the  coefficients  A. ..  ,   A..  _ ,   A-..  , 
A29    are   integrals  with  respect   to   time.      In  the  numerical   routine,    this 
integration   is   accomplished   by  a   Runge-Kutta  numerical    integration   scheme, 
contained  within  which   is   a  Newton-Raphson   iteration  to  generate  the  neces- 
sary changes    in  u,v,w  to   accomplish  the   integration. 

The  mathematical   basis    for  the  Newton-Raphson   iteration   scheme  to 
generate  corrections   for   u,   v,    w  is   the  following.      From  the  equations 

H        du   +  H        dv  +  H        dw  =   -   H 
uu  uv  uw  u 

(5.20)  h       du  +  H       dv  +  H        dw  =   -  H 

VU  VV  VW  V 
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we  find 


H  „  du  +  H   dv  +  H   dw  =  -  H 
wu       wv       ww         w 


du  = 


-  H        H 

u        uv 

H 

UW 

-  H        H 
V        w 

H 
vw 

~  \    *w 

\w 

(5.21) 


dv  = 


%u 

-  «u 

Huw 

H 

vu 

-  H 

V 

H 

vw 

™W1 

-ttw 

Hww 

D 


dw  = 


H 

uu 

H        -  H 

uv          u 

H 

vu 

H        -  H 

W              V 

H 

wu 

H        -  H 
wv          w 

where 


D  = 


H    H    H 

UU    UV    UW 


H    H    H 
vu   w   vw 


H  ,  H    H 

WU    WV    WW 


The  iteration  scheme  then  has  the  form 


(5.22) 


u„  =  u  ,  +  du 
n    n-1     n 


v_  =  v  ,  +  dv_ 
n    n-1     n 


w„  =  w„  n  +  dw 
n    n-1     n 


where  n  is  the   index  for  the  Newton-Raphson  iteration;   du,   dv,   dw  are 
those  found   in  equations    (5.21). 
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6.      Control    Variables  on  the  Boundary. 

In  section  4,    it  was   noted  that   for  a  path  to   be  admissible  the  con- 
trol variables   v,   w  had  to   satisfy  the   inequalities 

0   ±    w   £.   w 


(4.1) 


max 


Oi    vi    ' 


max 


Let  us  consider  the  situation  in  which  the  depth  w  assumes  the  max- 
imum depth  w    for  part  or  all  of  the  route  and  the  other  bounded  con- 

^    max     r 

trol  variable  v  remains  within  its  prescribed  bounds.  We  must  amend  the 
routine  for  determining  the  controls  as  follows.   We  take  the  maximum 

depth  w    of  the  type  of  submarine  being  considered  and  read  this  infor- 

r  max  ' r  ° 

mation   into   the  program  of  the  numerical    routine  which  calculates   the 

route.      If   the  Newton-Raphson   iterations   as   established    in    (5.22)   yield 

a  value   for  w  which  exceeds   w        ,    we   set  w  equal   to  w         and  generate 

max  max 

u  and  v,    using  the  following   iteration  scheme.      From  our   Newton-Raphson 
iteration  equations, 


(6.1) 


Huu  du  +  Huv  dv  =   "  Hu 


H       du  +  H       dv   =   -  H   , 

VU  W  V 


we  get  corrections  to  the  controls, 


-  H       H 
u       uv 


(6.2) 


du  = 


dv  = 


-  H       H 
V        w 


H  H 

uu       uv 


H         H 

vu       vv 


H  -  H 

uu  u 


H  -   H 

VU  V 


H  H 

uu  uv 

H  H 

vu  vv 
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Then  the  iterations  discussed  in  the  previous  section  take  the  form 

u  =  u  ,  +  du 
(6.3)      n    n"1 

Vn  =  vn-l  +  dvn  ■ 

Thus  we  may  modify  our  Newton-Raphson   iteration  scheme   for  the  con- 
trols when  the  minimum  value  of  H  occurs   for  w  on  the  boundary.      The 
modified  Newton-Raphson   iterations  generate   successive  corrections  to   u 
and  v  so   as   to  produce  an  admissible  path.      The  numerical  routine  is  pro- 
grammed  so   that   it   is   possible  for  the  depth  to   assume   its  maximum  for 
some  part  of  the  route  without   remaining   fixed  at  maximum  depth  after 
once  assuming   it.      In  section  7  we  will   see  an  optimum  path  which  has 
the  control  on  the  boundary  for  part  of  the  path;    the  submarine  travels 
at  maximum  depth  for   a  while,    then  comes  up. 

The  subroutine  for  calculating  u,   v  when  w  =  w         is  contained   in 

max, 

Appendix  I,  part  C.  The  above  features  can  be  seen  by  examining  either 

the  flow  chart  for  subroutine  BOUNDW  or  the  subroutine  itself  which  is 

given  in  part  G  of  Appendix  I. 

Modifications  similar  to  those  made  for  w  on  the  boundary  would  be 

made  if  v  =  v    for  some  part  of  the  route.   One  additional  change  re- 
max  r  to 

quired  in  this  case,  which  was  not  necessary  when  w  =  w   ,  is  that  dV 

be  set  equal  to  zero  in  the  numerical  routine  whenever  v  =  v    and 
^  max 

H   <  0.   V  =  v    and  H  >  0  imply  that  the  velocity  is  decreasing  or 
v  max      v  c  J  J  ° 

moving  away  from  the  bound  and  hence  the  numerical  routine  described  in 

equations  (5.22)  for  generating  corrections  to  the  control  variables 

would  be  used  whenever  this  is  the  case.   It  should  be  noted  that  setting 

(f\j   equal  to  zero  in  the  case  of  •>;  =  w    was  not  required  since  (fvr   does 

not  enter  into  the  equations  for  the  numerical  routine  in  section  5. 

When  v  =  v    and  H  <  0,  we  get  from  equations 
max      v 
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(6.4) 


H     du  +  H     dw  =  -  H 

uu              uw  u 

H     du  +  H     dw  =   -  H 

WU                    WW  W 


the  corrections  to  the  control  variables  u,   w  as 

-  H       H 
u       uwi 


du  = 


~  Hw     Hww 


(6.5) 


dw  = 


H  H 

uu  uw 


H  „,       H 

WU  WW 


H  -  H 

uu  u 


H  -  H 

WU  w 


^U        Huw 

WU  ww 


Using  the  du,   dw  found   in  equations    (6.5),   we  get   equations   for  correct- 
ing u,   w 


(6.6) 


u     =  u     ,    +  du 
n         n-1  n 


w„  =  w     .    +  dw 
n         n-1  n 


If  both  v  =  v    and  w  =  w    and  it  is  also  true  that  H  <  0  and 
max         max  v 

H  <  0,  we  generate  corrections  to  the  control  variable  u  by  using  the 


equation 

(6.7) 

from  which 
(6.8) 


H  du  =  -  H 
uu       u 


-  H 


du  = 


u 


H 


uu 


and  the  iterations  to  correct  u  take  the  form 

(6.9)        u_  =  u  .  +  du  . 

n    n-1  n 

Both  of  these  situations  were  encountered  in  generating  Path  IV 
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which  will  be  discussed  in  section  8.   The  results  given  there  will  indi- 
cate that  the  modifications  needed  when  v  =  v    or  v  =  v    and  w  =  w 

max        max         max 

do  not  affect  convergence  of  the  Newton- Raphson  iteration  schemes  in  the 
numerical  routine,  because  we  will  find  that  Path  IV  is  admissible. 

The  flow  chart  for  the  subroutine  BOUNDV,  which  is  used  when  v  = 
v  ,  can  be  seen  in  part  D  of  Appendix  I.  Part  B  of  Appendix  I  con- 
tains the  flow  chart  of  subroutine  VUW  which  takes  care  of  the  case  where 

v  =  v    and  w  =  w   .   The  result  of  setting  <fv   =  0  when  v  =  v    and 
max         max  max 

H  <  0  can  be  seen  by  examining  the  flow  chart  of  the  numerical  routine 
in  part  A  of  Appendix  I.   The  deck  listings  of  BOUNDV, VUW,  and  the  nu- 
merical routine  can  all  be  found  in  part  G  of  Appendix  I. 
7.   Paths. 

Our  computations  have  established  the  existence  of  three  extremals. 
The  three  paths  all  satisfy  the  Euler-Lagrange  equations  as  outlined  in 
section  4. 

These  three  paths  can  best  be  compared  by  listing  the  contrasting 
points  of  the  three.   The  areas  in  which  the  greatest  difference  appeared 
among  the  three  routes  are  the  probability  of  detection  p(T) ,  the  depth 
w,  and  the  constants  (h..,h  )  which  yield  admissibility. 

The  following  is  a  list  of  the  results  for  each  path  in  the  three 
areas  just  mentioned. 

Path  I 
p(T)  .15741  x  10~3 

w  thermocline  ~  200  feet 

(h1,h2)  (-.00019,  .00024) 

Path  II 
p(T)  .22725  x  10~3 
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w  w    up  to  68  3  feet 

max 

(h1,h2)  (-.00038,  .00048) 

Path  III 
p(T)  .23060  x  10"3 

w  525  to  725  feet 

(h1,h2)  (-.00037,  .00046) 

Path  I  gives  us  an  absolute  minimum,  i.e. ,  a  minimum  under  either 
weak  or  strong  variations.   Path  II  gives  us  a  relative  minimum  if  only 
weak  variations  are  allowed.   Path  III  is  an  extremal,  but  does  not  fur- 
nish a  relative  minimum  under  either  weak  or  strong  variations.   We 
called  Path  III  a  worsimax  path. 

If  any  additional  information  concerning  any  one  of  the  above  paths 
is  desired,  copies  of  the  three  paths  can  be  found  in  Appendix  II.   Given 
there  is  a  printout.,  of  each  path  with  the  coordinates  (x,y)  denoting  the 
submarine's  location,  the  control  variables  u,v,w,  and  the  probability 
of  detection  p(t)  for  each  time  step. 

Note  that  in  Path  II  w  =  w    for  nearly  all  of  the  route.   In  cal- 

max 

culating  this  path,  the  iteration  scheme  described  in  equations  (6.1), 
(6.2),  (6.3)  was  used.  The  fact  that  Path  II  converges  to  the  desired 
terminal  point  (x_,y  )  substantiates  the  results  given  in  section  6. 

Analysis  of  Paths.   Checks  on  the  conditions  as  outlined  in  section 
4  point  out  the  following  facts.   Path  I  satisfied  the  following  condi- 
tions: 

1.  admissibility  conditions, 

2.  Euler  equations, 

3.  Legendre  conditions,  and 

4.  Weierstrass  condition. 
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Path  II  satisfied 

1.  admissibility  conditions, 

2.  Euler  equations,  and 

3.  Legendre  conditions, 
but  not  the  Weierstrass  condition. 

Path  III  satisfied 

1.  admissibility  conditions,    and 

2.  Euler  equations, 

but  not  the  Legendre  conditions  nor  the  Weierstrass  condition. 

Hence,  by  the  criterion  established  in  section  4,  there  is  but  one 
minimum  path,  that  being  Path  I.   Note  that  the  check  of  the  envelope 
condition  was  not  included  in  this  investigation. 

The  above  results  point  out  an  important  fact  which  is  often  ignored 
or  overlooked;  the  generation  of  an  extremal,  by  no  means  guarantees  that 
you  have  the  desired  minimum.   This  emphasizes  the  need  for  a  check  on 
all  of  the  conditions  for  a  minimum  at  each  time  step  of  the  numerical 
integration  scheme  for  generating  the  path.   If  checking  after  each  time 
step  requires  excessive  computer  time,  the  checks  may  be  performed  at 
some  appropriate  periodic  intervals. 

It  can  be  noted  at  this  time,  that  if  we  restrict  ourselves  to  weak 
variations  as  defined  in  section  4,  both  Path  I  and  Path  II  are  extremals 
which  yield  relative  minima.   In  contrast  to  this,  analyzing  the  paths 
and  considering  strong  variations  yields  the  result  noted  above,  namely, 
that  Path  I  is  really  the  only  relative  minimum  among  the  three  paths. 

After  Path  II,  which  does  not  give  a  relative  minimum,  was  genera- 
ted, a  search,  as  outlined  in  the  discussion  of  the  Weierstrass  condition 
in  section  4,  was  used  to  determine  a  new  set  of  control  variables  u,v,w, 
which  would  satisfy  the  Weierstrass  condition.   The  resulting  paths 
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turned  out  to   converge  to   Path  I.      The  subroutine  search  was   used  for 
this  purpose  and  is  contained   in  Appendix  I. 

A  similar  search  over  the  grid  outlined   in  section  4  was   performed 
when  the  Legendre  conditions  were  not  met   in  the  case  of   Path  III.     The 
set  of  control  variables   which  gave  the  smallest  value  for  the  Hamilton- 
ian  in  this   search  were  then  used  to   continue  the  numerical   routine. 
Again  the   sequence  of  paths   produced  by  the  numerical    integration  con- 
verged to  Path   I. 

A  path  was    judged  admissible  if  it   came  within  one-fourth  mile  of 
the  desired  endpoint;    it  was   felt  that   further  accuracy  was   not  worth  the 
computing  time. 

To  test  the  convergence  of  the   numerical   routine,   on  a  few  paths   the 
routine  was   allowed  to   continue  until   no   further   improvement  occurred. 
In  each  of  the  three  paths   above,   duplication  occurred  with  accuracy  of 
at  least  one-tenth  of   a  nautical  mile.      Duplication  here  means  the  abil- 
ity of  the  numerical   routine   to   repeat    itself  after  once  converging  to 
the  desired  terminal   point. 

It  has  been  noted  that   a   spiral   pattern  of  convergence  about   the 
desired  terminal   point   is  present   in  the  computation  of  each  path.      It 
is  not  clear  why  this  occurs  but  the  following   is  offered   as  a  possible 
explanation.      In  the  computation,   we  take  H     =  H     =  H     =0  and  vary  u,v, 
w,x,y,   h,,h   ,   but  we  assume  that  x,y  have  the  values  they  assume  on  the 
path.      For  the  purposes  of  explanation,    let   us  examine  the  equation 
H     =0,    from  which  we  get   an  equation  of  the  form 

H  ,    dhn    +  H.    dh0   +  H,„<fu  +  H„„<Tv  +  H„  <fw  + 
uh,      1  uri2      2  uu  uv  uw 

H     <fx  +  H     <fy  =  0 

ux  uy  J 
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when  we  take  its  variations.  The  last  two  terms  in  this  expression  drop 
out  in  linear  problems  and  can  be  shown  to  be  negligible  if  T  is  small 
in  any  case.   They  introduce  considerable  complication  and  extra  computa- 
tion and  hence  were  discarded.   This  omission  may  be  the  reason  for  the 
spiral  pattern  of  convergence.   It  should  be  pointed  out  that  the  terms 
were  omitted  in  only  the  correction  routine.   If  the  sequence  of  paths 
converges,  there  is  no  related  error  in  the  path  to  which  they  con- 
verge. 

In  generating  Path  II,  a  subroutine  was  used  within  which  the  head- 
ing was  fixed  and  a  search  over  the  depth  and  the  speed  was  conducted  to 
determine  which  set  of  values  for  these  two  controls  gave  the  smallest 
value  for  the  Hamiltonian  TH.   These  values  were  then  used  in  the  numeri- 
cal routine  to  insure  a  start  in  the  proper  direction.   This  method  was 
put  into  use  when  it  was  found  that  poor  initial  choices  for  u,v,w  caused 
the  Newton- Raphson  routine  to  diverge  at  the  beginning  of  the  route.   This 
subroutine  SEARCH  can  be  seen  in  Appendix  I,  part  G/ 

To  insure  a  proper  start  in  computing  Path  II,  the  subroutine  WORSI 

was  used.   This  subroutine  is  the  same  as  subroutine  BOUNDW  described  in 

section  6  with  one  exception,  that  being  that  w  =  w    is  replaced  by 

max 

w  =  500  or  some  other  intermediate  value  for  the  depth.   This  subroutine 
then  fixes  w  at  500  and  computes  u  and  v  using  the  iteration  scheme  des- 
cribed in  equations  (6.3).   The  resulting  set  of  values  for  u  and  v  are 
then  combined  with  w  =  500  to  make  up  the  initial  guesses  for  the  numer- 
ical routine.   This  subroutine  is  also  given  in  Appendix  I,  part  G. 

Few  problems  were  encountered  in  the  generation  of  Path  I,  but  the 
introduction  of  conditions  to  cause  the  submarine  to  assume  its  maximum 

depth  w    or  velocity  v    aggravated  the  situation  and  introduced  dif- 
max  max 
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ficulties  to  make  the  subroutines  listed  above  necessary. 

It  should  be  realized  that  this  problem  contains  some  real  difficul- 
ties, if  approached  blindly.   With  the  proper  background  and  forethought, 
most  of  the  difficulties  can  be  anticipated  and  handled,  when  encountered, 
by  methods  such  as  those  described  above. 
8.   Corner. 

This  section  contains  a  discussion  of  a  case  in  which  the  path  gen- 
erated contains  a  corner. 

A  corner  appears  when  control  variables  u,v,w  which  minimize  the 
Hamiltonian  are  discontinuous  functions  of  t.  For  convenience,  the  nota- 
tion 


will  be  used  to  denote  a  set  u,v,w  of  control  variables.   The  conditions 
for  a  corner  are,  first,  that  there  exists  some  point  on  the  route,  call 
it  t,,  where  two  sets  of  control  variables  U  and  U  both  minimize  the 
Hamiltonian,  H,  and  second  that  one  set  of  control  variables,  say,  U  , 
gives  a  smaller  value  for  H  when  t  ^  t. ,  whereas  the  other  set  of  con- 
trols  U  yield  a  smaller  H  for  t  >  t,  .   These  two  conditions  can  be 
stated  as,  first 

(8.1)     HCU1)  =  H(U2)  =  min     t  =  t 


and  second 


(8.2) 


U 


H(UX)  <  H(U2)    for  t  ^  t 
H(U  )  >  H(U2)    for  t  >  t 


for  some  neighborhood  of  t  . 

At  a  corner  the  numerical  routine  for  the  corrections  should  be 
amended  by  adding  terms  to  handle  the  changes  due  to  the  variation  of 
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the  corner  time  t  . 

The  changes  necessary  are  outlined  in  the  earlier  report  "Optimum 
Submarine  Routing"  [l],  section  6.   Even  without  these  correction  terms, 
the  numerical  routine  generated  an  admissible  path  which  contained  the 
corner,  namely  Path  IV. 

To  construct  a  mathematical  model  in  which  a  corner  could  be  ex- 
pected the  functions  described  in  section  2  were  changed  in  the  following 
manner.  We  replaced 

.  1  w 


1  + 


w0 


1  + 


in  the  equation  for  f     by 


w 


W0 


1 


w 
1  +  

2 
wl 

with  w,  equal  to  five  hundred  feet.   The  constants  in  f,,  f_,  the  equa- 
tions representing  the  passive  and  active  defenses  respectively,  were 
chosen  in  such  a  way  that  a  corner  could  be  anticipated. 

The  model  used  was  one  in  which  the  passive  defense  was  dominant 
for  the  beginning  of  the  route,  the  two  would  become  equal  at  approxi- 
mately the  middle  of  the  route,  and  the  active  search  would  then  dominate 
for  the  latter  part  of  the  route.   These  facts  are  apparent  when  one  looks 
at  Path  IV  in  Appendix  II. 

When  we  analyze  Path  IV,  we  see  that  the  submarine  proceeds  at  ap- 
proximately thermocline  depth  for  the  first  part  of  the  route  and  then 

changes  to  w  =  w    for  the  remainder  of  the  route.   It  can  also  be  noted 
°  max 

that  the  speed,  v,  was  considerably  less  than  v  =  v    until  the  corner 

max 
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was  encountered,    at  which  time  v  became  equal    to  v        .     Whenever  the  active 

max 

search  completely  dominates,  the  submarine  goes  as  deep  and  as  fast  as 
possible. 

A  check  on  the  Hamiltonian,  H,  after  each  time  step  shows  that  H  is 
constant,  within  the  accuracy  of  the  routine,  for  the  time  steps  before 
we  reach  the  corner,  but  is  not  constant  as  we  proceed  beyond  the  corner. 
It  is  not  clear  why  H  does  not  remain  constant  throughout  the  route,  but 
a  possible  explanation  is  that  the  corner  was  effectively  passed  before 
it  was  found,  i.e.,  the  numerical  routine  failed  to  detect  the  corner 
when  the  conditions  for  a  corner  were  in  fact  present.  The  failure  to 
find  the  corner  immediately  is  a  result  of  the  grid  size  used  in  the  sub- 
routine SEARCH  to  compare  the  Hamiltonian  for  controls  U  generated  by  our 
numerical  routine  with  the  controls  U  used  in  the  search.   As  noted  be- 
fore, this  grid  size  was  decided  upon  when  a  finer  grid  was  found  to  re- 
quire excessive  computer  time.   Considering  that  it  takes  a  while  for  the 
routine  to  find  the  corner  and  it  takes  the  routine  a  certain  amount  of 
time  to  settle  down  after  the  corner  is  found,  the  fact  that  admissibil- 
ity was  accomplished  was  thought  sufficient  to  justify  omission  of  the 
corner  correction  terms. 

Notice  that  in  Path  IV  v  =  v    with  w<w    for  a  part  of  the  route 

max         max 

before  the  corner  and  then  v  =  v    and  w  =  w    for  the  portion  of  the 

max         max         c 

route  that  comes  after  the  corner.   Path  IV  was  generated  using  the  iter- 
ation schemes  described  in  equations  (6.4),  (6.5),  (6.6)  when  v  =  v 

max 

and  w  <  w    and  by  using  (6.7),  (6.8),  (6.9)  when  v  =  v    and  w  =  w 

max  max         max 

The  admissibility  of  Path  IV  confirms  the  results  stated  in  section  6  for 
v  on  the  boundary  and  when  both  v  and  w  are  on  the  bound  of  their  allow- 
able values. 
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9.      Observations. 

This   section  contains  observations  which  may  be  helpful   to   a  person 
wishing  to   continue  the  study  of   the   submarine  routing  problem. 

In  the  Newton-Raphson  iterations  which  occur  it  may  be  possible  to 
improve  both  convergence  and  accuracy  as  follows.  For  example,  in  gen- 
erating the  control  variables   let  us  make  up  an  error   function 

2  2  2 

e.    H       +  e0  H       +  e   'H     . 
1      u  2     v  3     w 

If  each  successive  iteration  does  not  diminish  this  error  function,  then 
we  should  diminish  the  preceeding  corrections  by  a  factor  of  say,  two, 
or  five.   The  reason  is  that  the  Newton-Raphson  iteration  moves  the  var- 
iables in  the  right  direction  but  may  overshoot  if  the  linear  terms  are 
not  dominant.   The  iteration  would  be  terminated  when  the  above  error 
function  was  less  than  some  preassigned  value.  The  incorporation  of  such 
routines  might  well  improve  convergence,  save  computing  time,  and  improve 
accuracy.   The  convergence  criterion  above  is  derived  from  the  fact  that 
satisfaction  of  the  Euler  equations  implies  H  ,H  ,  H  are  all  equal  to 

^  U   V    w 

zero.   Similar  conditions  could  be  established  for  the  other  Newton- 
Raphson  iterations  for  generating  corrections  to  the  parameters  h,,tu. 
In  the  iterations  to  correct  these  the  established  criterion  would  be 
based  upon  admissibility  of  the  path.   By  letting  (x,y)  represent  the 
endpoint  of  the  path  generated  and  (xT,yT>  be  the  desired  terminal  point, 
we  could  write  the  condition  as 

(x  -  xT)   +  (y  -  yT)  <  i    . 
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10.   Summary. 

In  this  paper  the  submarine  routing  problem  was  studied.   Functions 
were  chosen  that  seemed  to  be  typical  of  the  functions  representing  the 
detection  devices,  both  passive  and  active.   Information  that  could  be 
arrived  at  only  through  the  use  of  empirical  data  such  as  the  sea  state, 
for  example,  was  not  made  a  part  of  these  functions. 

Using  this  mathematical  model  and  determining  the  path  from  (xQ,yQ) 
to  (xT,yT) ,  in  a  fixed  time  T,  which  minimizes  the  probability  of  detec- 
tion, p(T) ,  resulted  in  the  generation  of  three  extremals.   Examination 
of  the  paths  using  established  criterions  for  a  relative  minimum  lead  to 
the  following  results:  one  path  yielded  the  desired  minimum,  a  second 
satisfied  all  conditions  except  the  Weierstrass  condition  and  the  third 
path  was  just  an  extremal,  satisfying  neither  the  Legendre  nor  the  Weier- 
strass conditions. 

Situations  were  encountered  in  which  the  speed,  v,  or  the  depth,  w, 
or  both  v  and  w  were  on  the  boundary  of  allowable  controls.  It  was  found 
that  if  the  control  on  the  bound  was  set  equal  to  the  boundary  value  and 
corrections  generated  for  the  remaining  control  variables,  admissibility 
was  accomplished  just  as  when  all  controls  were  interior  to  their  allow- 
able ranges. 

With  a  change  in  the  original  model  for  the  active  defense,  condi- 
tions conducive  to  a  corner  were  established.   The  numerical  routine 
then  generated  a  path  which  had  a  corner  and  was  admissible.   Admissi- 
bility here  was  accomplished  without  the  use  of  corrections  for  the 
corner,  and  for  this  reason  the  corrections  were  not  made  a  part  of  the 
numerical  routine. 

The  numerical  routine  as  described  in  section  5  was  programmed  in 
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Fortran  1960  for  a  CDC  1604  computer.   Both  a  flowchart  and  the  program 
for  this  routine  can  be  found  in  Appendix  I,  the  flow  chart  in  part  A 
and  the  deck  listing  for  the  program  in  part  G. 
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Appendix  I 

This  appendix  contains  the  f}.ow  charts  and  the 
deck  listings  for  the  program  of  the  numerical  routine 
and  the  subroutines  that  were  used  to  generate  the 
paths  described  in  the  text. 
A,  Flow  chart  for  the  numerical  routine. 

(START?) 

vk 


Dimension  YVARS,YC,DY,C,AK,TAU,X,Y,Z,UT,VT,WT 


sy 

/^Read  XT,T,TEETA,A1,B1,C1,D1,S1,C2,A 

D2,X2,Y2,H1,K2,FaC,FMUL,W0,LMAX, 
VJ^VJMAX.V^-:.      _1_  / 


-JS£__„ 


JQ  =  THEIA  /  57.29, 
COS  =  COSF(Q) 
SIM  =  SIKF(Q) 

X2  =  SI  COS  +  XT 

12  =  SI  SIM 
VAV  =  XT  /  T 


v. . 

Print  XT,T,TSE!£A,Al,Blt 01,01,31,02/ 
J  D2,X2,Y2,H16H2,FAC,F^IDL,W0,LMAX, 
VW1,MMAX,VMAX 


!  C(l)  =  0.0   .  in  =  T  /  XSTEP  j 
:  C(2)  =  0.0      XN1  =  Ml 
C(3)  =  0.0    STEP  =  T  /  XN1 
C(4)  =  0.0     N2  =  Nl  +  1 

!  XSTEP  =  T  /  FAC I 




U  =  -  .2| 
V  =  VAV  ! 

w  =  wo 
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IC(J)  =  YVARS(J)  +  C(I)  AK(I  -  1,J),J  =  1,9 


_i2_ 


Compute  DBTER,gUMl,FUM2,FUH3 


> 
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1:1 


> 


uzn,RO  =  u; 

VZiRO  =  V ! 
WZERO  =  W 


N/ 


jCompute  FlT, tt2, FX.FY , Sll , S12 , S21 , S22]e 


2-C3)  =  COSU  S21  -  V  SINU  Sif 

BT (4)  =  COSU  S22  -  V  SINU  S12 

|DY(5)  =  SINU  S21  +  V  COSU  Sll 

DY(6)  =  SINU  S22  +  V  COSU  S12J 


DY(3)  = 
DY(4)  = 
DY(5)  = 
DY(6)  = 


-  V  SINU  Sll 

-  V  SINU  S12 

V  COSU  Sll 

V  COSU  S12 


pl(7)  =  FI1  +  FI2 
M  DY(3)  =  V  COSU 
LDY(9)  =  V_SINU__ 


~>; 


-±&L 


!AK(ltJ)  =  STiSP  DY(J),J  a  1,9 


|YVABS(J)  -  YVARS(J)  +  (AK(I,J)  +  2  AK(2,J)  +j 
12  AKQ.J)  +  AK(4,J))   /  6, J  =  1.9 __J 
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TVAR  =  TVAR  +  STEP 

TAU(K)  =  TVAR 

X(K)  =  YVARS(8) 

Y(K)  =  YVARS(9) 

Z(K)  ■  YVARS(?) 

UT(K)  =  U 

VT(K)  =  V 

V;T(K)  a  W 


"DET  =  YVARS(3)  XVARS(6)  -  YVARS(5)  YVARS(4) 


IXNUMl  =  (XT  -  YVARS(8))  YVARS(6)  +  YVARS(9)  YVARS(4-) 
JXNUK2  =  -  YVARS(9)  YVaRS(3 )  -  (XT  -  YVARS (8))  YVARS(5) 


.^L. 


HI  =  HI  +  FMUL  XNUK1  /  D2T 
H2  «  H2  +  FMU  ,  XUDM2  /  PET 


(  Print  K2  ,H1 ,H2~Tx(H2) , Y(N2) t Z(K2~7) 


(Print  TaU(I)>X(I).Y(I)>DT(I),VT(I),VIT(I)"72(I).I  =  1.N2  ) 


IMAX^>~ ^L  -  L  +  i 


PROS  =  1.0  -  SXPF(-  0.0001  2(N2)) 


(Print  FROB 


STOP  I 
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3.  Flow  chart  for  subroutine  VUW. 


0- 


I  START  1 

51=3 


DU 

=  FUM1  /  DETER 

DV 

=  FUM2  /  DETER 

DW 

=  FUK3  /  DETER 

.^L 


U  =  D  +  DU 
V  =  V  +  DV 
W  =  W  +  DW 


Comnute  FUM1,FUM2J< 

FUM3,  DETER 


\k 


RETURN  fc-j 


->:i  =  i  +  i; 


:\ 


JJ   Call 
-.     'iBOUNDV 


'  <HV:0: 


-^|K  =  1 


^k: 


«- 


[Compute  HU.HUU 


<*) 
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C*  Flow  chart  for  subroutine  BOUNDW. 


START] 


Jl. 


W  =  WAX 


\/ 


Compute  HU,HV,KUU,HVU,HW 


0- 


y/- 


D2N  =  HUU  HW  -  HVU  HVU 
GUM  =  -  HU  HW  +  HV  HVU 
GUM2  =  -  HV  HUU  +  HU  HVU 


i 


1  =  1 


|DU  =  GUM1  /DEN 
[DV  =  GUM2-  /  PEN 


-V- 


U  =  U  +  DU 
V  =  V  +  DV 


V 

| Compute  HU, HV , HUU , HVU, HVV] 

i 


[DEN  =  KUU  HW  -  HVU  HVU] 
iGUMl  =  -  HU  HW  +  HV  HVU| 
JGUM2  =  -  HV  HUU  +  HU  HVUj 


\1/ 

[rEtuMI 


-^1  =  1  +  11 


0 
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D.  Flow  chart  for  subroutine  BOUNDV. 


!  START]' 


(V  =  VH.4X1 


„\z 


Compute  EU,KW,HUU,HWU,HWW 


[BEET  ■  HDD  HWW  -  HWU  HWU 

JGDM1  =  -  ED  HWW  +  HW  HWU 


IGUK2  =  -  HW  HUU  +  HU  HWU 


Ul 
Hi 


\/ 


& 


1  =  1 


DD 

DW 


GUM1  /  DEM i 
GDM2  /  DSNj 


^L 


U  =  U  +  DD 
W  =  W  +  DW 


V 


(Compute  HU,HW;hUU,HWU7BWW1 


v 


]dsn  =Thdd  kww  -  ewu  hwu 

j'GUMl  =  -  HU  HWW  +  hw  hwu 
JGUM2  =  -  ffi^J  HUU  +  HU  HWU 
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E.  Flow  chart  for  subroutine  SEARCH. 


jSTART] 


Dimension  VIKT,WINTtH_ 


V  . 


HAM  ■  XLAM  V  COSU  +  XMU  V  SINU  +  FI1  +  FI2 

T 


\/ 


VDJT(O)  =  0.0! 
WIHT(O)  =  0>0l 


I  VINT(M)  =  VINT(M  -  1)  +  3.0 


d> 


^/.. 


N  =  1 


[vmgrCK)  =  MI(N  -  1)  +  200.0] 


jAl 


JH(M.N)  ■  XLAM  VDiT(M)  COSU  +  XMU  VIKT(M)  SINU  +  FI1  +  ml 

I 


< 


^[N  =  N  +  1! 


/■£. 


^M  =  M  +  1 


HAMSRCH  =  H(l,l) 
YSRCH  =  VIKT(i) 
VJSRCH  =  WINT(l) 


XD 


O 


48 


k9 


F.  Flow  chart  for  subroutine  WORSI. 


!W  =  500. I 


^k- 


Compute  KUtHV.HUU,HVU,HW 


J^_ 


DEN  =  HUU  HW  -  HVU  HVU 
GUM1  =  -  HU  HW  +  HV  HVU 
GUM2  =  -  HV  HUU  +  HU  HVU 


1 


DU  =  GUI41  /  DEN 
DV  =  GUI-I2  /  DEN 


\/ 


!U  =  U  +  DU 

!v  =  v  +  dv 


\!/ 


Compute  HU.HVtHUU.HVU.HWl 


^L 


DEN  =  HUU  HW  -  HVU  HVU 
GUM1  =  -  HU  HW  +  HV  HVU 
GUM2  =  -  HV  HUU  +  HU  HVU 


Z. 


-£1 


I  +  1 


1? 


G„.     Fortran  Program:     Printout  of  deck  of  punch  cards 
submitted  to  computer. 

The  main  numerical  routine  and  various  subroutines 
are  given  here.     The  equations  in  the  numerical  routine 
and  the  first  four  subroutines  as  listed  here  are  the 
ones  which  gave  a  corner.     The  equations  in  subroutine 
WORSI  are  the  ones  used  in  generating  the  three  extremals. 


PROGRAM    SU3R0UTE 
C       YVARS( I )=LAM3    YVARS(2)=MU3       YVARS(3)=A11       YVARS(4)=A12       YVARS(5)=A21 
C       YVARS(6)=A22       YVARS(7)=       Z       YVARS(8)=       X       YVARS ( 9 ) =       Y 
C       XLAM=H1+LAM3       XMU=H2+MU3 

DIMENSION    YVARSC9) >YC(9)>DY(9),C(4)>AK(4»9) »TAU(400) .X(900) ♦ 
+  Y(90     ),Z(400)»VT(400) tUT(400)»WT(400) 

READ    1,XT»T,THETA*A1»B1,C1»D1,SI»C2»D2*X2»Y2*H1»H2,FAC»FMUL» 
+W0 , LMAX  » W 1 , WM AX  >  VM A X 

1  FORMAT     (:r5.0»2F4.0*F2.0*F3.2»F?.l,F5.4,F4.0»F7.6»F5.4»F6.1»F4.1» 
+  F6o4,F6.0>F4.0,F4.2»F4.0»I2/F4.0,F5.0,F3.0) 

Q=THETA/57.2957  7951 

CC5=CCSF(Q) 

SlN=SfNF(Q) 

X2=S1*C0S+XT 

Y2=S1*SIN 

VAV=XT/T 

PRINT  2 

2  FORMAT ( 1H02HXT4X1HT3X5HTHETA2X2HA12X2HB12X2HC12X2HD12X2HS12X4HWMAX 
+1X4HVMAX) 

PRINT  3>XT,T»THETA,A1,B1,C1,D1,S1»WMAX,VMAX 

3  FORMAT  (3F5«0»F5.1»F4.2,F4.1,F5.4,F5.0,F5.0»F3.0/) 
PRINT  4 

4  -ORMAT(3H  C24X2HD24X2HX24X2HY2 5X2HH17X2HH23X3HFAC2X4HLMAX 
+3X4HFMUL2X2HW05X2HW1 ) 

PRINT  5>C2»D2»X2>Y2»H1»H2»FAC,LMAX,FMUL,W0»W1 

5  FORMAT  (F6.5»F5.1»F6.1»F4.1,2F9.5»F6.0,I3»F8.2»F6.0»F6.0/) 
C   DEFINE  RUNGE-KUTTA  CONSTANTS. 

C(1)=0.0 
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C(2)=0.5 

C(3)=0.5 

C(4)=1.0 

X( 15=0.0 

Y( 15=0.0 

TAU(1)=0.0 
C   COMPUTE  TIME  STEP  LENGTH. 

XSTEP=T/FAC 

N1=T/XSTEP 

XNl=Nl 

STEP=T/XN1 

N2=N1+1 
C   GUESS  INITIAL  VALUES  FOR  U,  V,  AND  W  BEFORE  FIRST  ITERATION 

U  =  -.9 

V  =  VAV 

W  =  WO 
C   ENTER  LOOP  FOR  GENERATING  CORRECTIONS  TO  H1»H2. 

DO  14  L=1,LMAX 

TVAR=0.0 
C   GUESS  INITIAL  VALUES  FOR  U»V»W,  AFTER  FIRST  ITERATION 

IF  (L-l)  198,199.19s 

198  U=UZERO 
V=V2ERO 
W=W2ERO 

199  UT(  D=U 
VT(1)=V 
WT(D=W 

DO  6  1=1,9 
6  YVARS{ I )=0o0 
C   ENTER  LOOP  FOR  COMPUTING  THE  PATH  FOR  EACH  TIME  STEP. 
DO  11  K=2,N2 
IFU-3)  901,902,902 

90  2  CALL  SEARCH(A1,31,C1,D1,FI1,FI2»XLAM,XMU,C0SU,U,V»W»SINU,W0» 
+G1,G2,W1,C0SUMQ) 
C   ENTER  LOOP  FOR  RUNGE-KUTTA  INTEGRATION. 
901  DO  88  1=1,4 
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DO  7  J=l,9 

YC( J)=YVARS( J)+C( I )*AK( I-1,J) 

XLAM=H1+YC( 1) 

XMU=H2+YC(2) 

XD=(XT-YC(8)  )*C0S+S1-YC(9)*SIN 

FIl=EXPF(-(XD*LOGF(2. )/1000.  >  ) 

FI2=C2+D2*EXPF(-( ( YC ( 8 J-X2 )**2 

COSU=COSF(U) 

SINU=SINF(U) 

COSUMQ=COSF(U-Q) 

SINUMQ=SINF(U-Q) 

G1=C1+D1* ( W-WO) *( W-WO ) 

G2  =  1 .+D 1* ( W-WO ) * ( W-WO ) 

G3=( .01)*FI1*( A1+B1*V*V) 

G4=(-G1*(2«*D1*(W-W0) )+G2*(2«* 

G5=(C1+D1*( W-WO) *( W-WO) )/( l.+D 

G7=1.+W*W/(W1*W1) 

Pl=( ( l.+Dl*(W-WO)*(W-WO) )*2.*D 

P2=(-(8.*Dl*(W-WO)*(W-WO> )*D1) 

P3=  (-(Cl+Dl*(W-WO)*(W-WO) )*2* 

P4=  ( Gl*8. *D1*( W-WO >*D1*( W-WO) 

P9=FI2*(-2./(Wl*Wl*G7*G7)+8.*W 

HV=XLAM*COSU+XMU*SlNU+( .01 J*FI 

HU=-XLAM*V*SINU+XMU*V*COSU+G3* 

HVV=( .01)*2.*Bl*FIl*G5*(l.-.25 

HVU=-XLAM*SINU+XMU*COSU+G5*.0l 

HUU=-XLAM*V*CGSU-XMU*V*S I NU+G3 

-  -G3*< l.-.2  5*(COSUMQ) )*G4+FI2 

-  U=G3*( .25*<SINUMQ) )*G4 
HWV=G4*( •01)*FI1*( 2«*B1*V)*(1. 
HWW=   G3*( ( l.-.25*COSUMQ)*(Pl+ 
P5=HVV*HWW-HWV*HWV 

P6= ( -HV ) *HWW+HW*HWV 
P7  =  HW*(-HW)+HWV*HV 
DETER=HUU*P5-HVU*(HVU*HWW-HWU* 
FUM1=(-HU)*P5-HVU*P6+HWU*(  (-HV 


+(YC(9)-Y2)**2)*LOGF(2.)/100000.) 


Dl*(W-WO) ) )/(G2*G2) 
1*( W-WO)* (W-WO) ) 

1+4.*D1*(W-W0)*D1*(W-W0) )/(G2*G2) 

/(G2*G2) 

*D1-4.*D1* (W-WO) *D1*( W-WO) )/(G2*G2) 

)/(G2*G2*G2) 

*W/(W1*W1*W1*W1*G7*G7*G7) ) 

1*(2.*B1*V)*G5*( l.-.2  5*COSUMQ) 

G5*( •25*SINUMQ) 

*COSUMQ) 

*FI1*2.*B1*V*( ,25*SINUMQ) 

*G5*( .25*COSUMQ) 

*(-2.*(W/Wl)*Cl./Wl) )/(G7*G7) 

-.25*(COSUMQ) ) 
P2+P3+P4) )+P9 


HWV  ) +HWU*  (  HVU*HWV-HWU*HVV  ) 
)*HWV+HW*HVV) 
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FUM2=HUU^P6+HU*(HVU-HWW-HWU*HWV)+HWU*(HVU*(-HW)+HWU*HV) 
FUM3  =  HUU*P7-HVU*(HVU*(-HW)+HWU*HVJ-HU*(HVU*HWV-HWU*HW) 

rAi i  SUBROUTINE  FOR  CORRECTIONS  TO  U»V»W. 

CALL  VUW(  DETER,  FUM1 .FUM2.FUM3,XLAM.XMU,FI1,FI2.A1.B1.C1.D1.W0. 

+W1»WMAX,VMAX»U.V»W,Q) 

FIV=COEFF  OF  FI1      F2V=COEFF  OF  FI2 

IF(K-2)  195, 297,19.5 
297  IF  (1-1)  195,196,195 
196  UZERO=U 

VZERO=V 

WZERO=W 
19  5  XD=(XT-YC(8) )*COS+Sl-YC ( 9 )*SIN 

FI  l  =  EXPF(-(XD*LOGF(2..  )/1000.  )  ) 

Fi2=C2+D2*EXPF(-((YC(8)-X2)**2+(YC(9)-Y2)**2)*LOGF(2.)/100OOO.) 

G1=C1+D1*( W-WO) *( W-WO) 

G2=l .+D1* ( W-WO ) * ( W-WO ) 

G3=( .01)*FI1*( A1+B1*V*V) 

G4=(-G1*(2.*D1*(W-W0l )+G2*(2.*Dl*(W-WO) ) )/(G2*G2) 

G5= ( C 1 +D1 # ( W-WO ) * ( W-WO ) ) / ( 1 . +D 1* ( W-WO ) * ( W-WO ) ) 

G7=1.+(W/W1)**2 

COSU=COSF(U) 

SINU=SINF(U) 

SINUMQ=SINF(U-Q) 

Pl=  (U~+D1*  (W-WO)  *<  W-WO)  )*2.*D1+4.*D1*<  W-WO  )*Dl*<  W-WO)  )/<G2*G2) 

P2=(-( 8.*Dl*(W-WO)*(W-WO) )*D1)/(G2*G2> 

P3=  (r-(Cl+Dl*(W-WO)*(W-WO) )*2.*Dl-4.*Dl*(W-WO)*Dl*(W-WO) )/(G2*G2 

P4=" ( G1*8.*D1*( W-WO )*D1*C W-WO) ) /(G2*G2*G2) 

P9=FI2*(-2./(Wl*Wl*G7*G7)+8.*W*W/{Wl*Wl*Wl*Wl*G7*G7*G7) ) 

HV=XLAM*COSU+XMU*SINU+(.01)*FI1*(2.*B1*V)*G5*(1.-.25*COSUMQ) 

HU=-XLAM*V*SINU+XMU*V*COSU+G3*G5*( .2  5*SINUMQ) 
HVV=(.01)*2.*B1*FI 1*G5*( 1 .-.2  5*COSUMQ ) 

HVU=-XLAM*SINU+XMU*COSU+G5*.0l*FIl*2.*Bl*V*( .2  5*SINUMQ) 
HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .2  5*COSUMQ) 
HW=G3*( l.-.2  5*(COSUMQ> ) *G4+FI 2* ( "2 .* ( W/Wl ) * ( 1./W1 ) ) /(G7*G7) 
HWU=G3*( . 25*(SINUMQ) )*G4 


I 
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HWV=G4*( .01 )*FI1*( 
HWW«  G3*<t!**-.25* 
F1V=(A1+81*V*V)*( . 
F2V=1./(1.+(W*W/(W 
FY*F1V*FI1*SIN*L0G 

+(YC(9)-Y2)**2)*LOG 

+  )  ) 
FX=F1V*FI1*SIN*L0G 

+(YC(9)-Y2)**2)*LOG 

+  )  ) 
DENOM=HVU* ( HVU*HWW 

+HWV*HVU ) 
511= ( <HW*HWW-HWV* 
Si  2= ( ( HVV*HWW-HWV* 
S21=( (HVU*HWW-HWU* 
S22=( t HVU*HWW-HWU* 
DY(1)=-FX 
DY(2)=-FY 

IF(V-VMAX)118»119» 
DY(3)=COSU*S21-V*S 
DY(4)=C0SU*S22-V*S 
DY(5)=SINU*S21+V*C 
DY(6)=SINU*S22+V*C 
,60  TO  120 

' V=VMAX 
DY(3)=         -V*S 
DY(4)=         -V*S 
DY(5)=  V*C 

DY(6)=  V*C 


2.*B1*V )*(!.-. 2 5*(0OSUMQ) ) 

C0SUMQ)*(P1+P2+P3+P4) )+P9 

01 )*(!.-. 2 5*C0SUMQ)*G1/G2 

1*W1  )  ) ) 

F( 2. >/1000.+F2V*(-D2*EXPF(-{ C YC ( 8 ) -X2 ) **2+ 

F( 2. )/ 100000.  )*(2.*(YC(9)-Y2)*LOGF(2.)/100OO0. 

F(2. )/l000.+F2V*(-D2*EXPF(-( ( YC ( 8 ) -X2 ) **2+ 

F( 2.  J/100000.  )*(2.*(YC(8)-X2)*LOGF(2.)/10O000. 

-HWU*HWV ) -HUU* ( HVV*HWW-HWV*HWV ) +HWU* ( HVV*HWU- 


120 


3 
88 


HWV ) * ( -V*S I NU )  +  < HVU*HWW-HWV*HWU ) * ( -COSU ) 
HWV ) * ( V*C0SU ) + ( H VU*HWW-HWV*HWU ) * ( -S I NU ) ) 
HWV ) * ( V*S I NU ) + ( HUU*HWW-HWU*HWU ) *COSU ) /DE 
HWV ) *  U -V ) *C0SU )  +  ( HUU*HWW-HWU*HWU ) *S I NU ) 


119 

INU*S11 

INU*S12 

0SU*S11 

0SU*S12 


)/DENOM 
/DENOM 
NOM 
/DENOM 


DY(7>=FII*F1V+FI2* 

DY(8)=V*C0SU 

DY(9)=V*SINU 

DO  8  J=i»9 

AK( I»J)=STEP*DY(J) 

CONTINUE 

DO  9  J=l»9 


INU*S11 
INU*S12 
0SU*S11 
0SU*S12 
F2V 
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9  YVARS(J 
TVAR=TV 
TAU(K)= 
X(K)=YV 
Y(K)=YV 
Z(K)=YV 
UT(K)=U 
VT(K)=V 
WT(K)=W 

11  CONTINU 
DET=YVA 
XNUM1=( 
XNUM2-- 
H1=H1+ 
H2=H2+ 

.  PRINT  1 

12  FORMAT 
PRINT  1 

13  FORMAT 

18  PRINT  1 

19  FORMAT 
PRINT  2 

20  FORMAT 

14  CONTINU 
PROB=l. 
PRINT  2 

21  FORMAT 
STOP 
END 


)=YVARS( J)+(AK( 1,J)+2.*AK(2>J)+2,*AK( 3 ♦ J ) +AK { 4, J ) )/6. 

AR+STEP 

TVAR 

ARS(8) 

ARS(9) 

ARS(7) 


E 

RS 
XT 
YV 


2 

(2 

3, 

(I 

9 

(1 

0, 

(F 

E 

0- 

1* 

(1 


(  3  )  *-YVARS  (  6  )  -YVARS  (  5  )  * YVARS  (  4  ) 
-YVARS ( 8 ) ) * YVARS ( 6 ) +YVARS ( 9 ) *YVARS ( 4 ) 
ARS ( 9 ) *YVARS ( 3 )- ( XT-YVARS ( 8 ) ) *YVARS ( 5 ) 
FMUL*XNUM1/DET 

FMUL*XNUM2/DET 

X2HN25X2HH17X2HH2  5X5HX(N2 ) 4X5HY ( N2 ) 6X5HZ ( N2 ) ) 
N2,H1»H2,X(N2 ) ,Y(N2 ) ,Z(N2) 
4,2F9.5,2F9.1,E13.5/) 


H0  3X3HTAU6X1HX7X1HY7X1HU5X1HV5X1HW10X1HZ/) 
(TAU( I ) »X( I ) »Y( I ) ,UT( I ) ,VT( I ) »WT( I) »Z( I) »I 
8.2»2F8«1»3F6«1»E13«5) 

EXPF(-0.0001*Z(N2 ) ) 
PROB 
H02X5HPROB=E13.5/ ) 


ltN2) 


SUBROUTINE  VUW( DETER *FUM1 ♦ FUM2 ♦ FUM3 »XLAM»XMU »FI 1 t FI 2 .Al , Bl »C1 »01 ♦ 
+WO»Wl»WMAX»VMAX»U»VtW»Q) 
DO  5  1=1,7 
DU=FUM1/DETER 


DV=FUM2/DETER 

DW=FUM3/0ETER 

U=U+DU 

V=V+DV 

W=W+DW 

IF(WMAX  -W)  390,390,392 
390  IF(VMAX-V)  1000,1000,1001 
392  IF(VMAX-V)  3000,3000,3001 
300  0  CALL  30UNDV(U,V,W,XLAM,XMU,FU,A1,B1,C1,D1,W0,FI2,Q,VMAX»HV»W1) 

IF(HV-C.O)  4000,4000,3001 
4000  RETURN 

1000  DO  2000  1=1,5 
V=VMAX 
W=WMAX 
CCSU=C0SF(U) 
STNU=SINF(U) 
COSUKQ=COSF(U-0) 
SINUMQ=SINF(U-Q) 
G1=C1+D1*{W-W0)*(W-W0) 
G2=l«+Dl*(W-WO)*(W-Wb) 
G3=C..01)*FI1*(A1+B1*V*V") 

G4=(-G1*(2«*D1*<W-W0)  )+G2*(2.*Dl*(W-WO)  )  )/(G2*G2) 
G5=(C1+D1*(W-W0)*(W-W0)  5/  (  1  .+D1*  ( W-WO)*  ( W-WO  )  ) 
HU=-XLAM*V*SINU+XMU*V*C0SU+G3*G5*< .25*SINUMQ) 
HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*(  .25*C0SUMQ) 
DU=-HU/HUU 

U=U+DU 
2000  CONTINUE 
RETURN 

1001  CALL    B0UNDW(U,V,W,XLAM,XMU,FU,A1,B1,C1,D1,W0,FI2»W1,0,WMAX,HW) 
IFIHW-O.O)     5000,5000,3001 

5000    RETURN 

3001    C0SU=C0SF(U) 

SINU=SINF(U) 

C0SUiviG  =  C0SF(U-Q) 

SINUMQ=SINF(U-Q) 
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G1=C1+D1*(W-W0)*(W-W0).  •  • 

G2 = 1. +D 1*( W-WO )*( W-WO) 
G3=(.01)*FI1*(A1+B1*V*V) 

G4=(-G1*( 2.*D1*(W-W0) )+G2*(2.*Dl*(W-WO) ) )/(G2*G2) 

G5=(C1+D1*(W-W0)*( W-WO) )/( l.+Dl*( W-WO)* (W-WO) ) 

G7=1.+W*W/(W1*W1) 

Pl=(  (l.+Dl*( W-WO)* (W-WO)  )*2.*Dl+4.*Dl*(W-WO)*Dl*(W-WO)  )/(G2*G2) 

P2=(-(  8.*D1*( W-WO >*( W-WO)  )*D1  )  /(G2*G2  ) 

P3  =  (-(Cl+Dl* (W-WO)* (W-WO) ) *2. *D1-4.*D1*( W-WO )*D1*( W-WO) )/(G2*G2) 

P4  =  (G1*8.*D1*(W-W0)*D1*( W-WO) ) /(G2*G2*G2) 

P9=FI2*(-2./(Wl*Wl*G7*G7)+8.*W*W/(Wl*Wl*Wl*Wl*G7*G7*G7) ) 

HV=XLAM*COSU+XMU*SINU+( .01 ) *FI 1* ( 2 .*B1*V ) *G5*( 1 .-.25*COSUMO ) 

HU  =  -XLAM*V*SINU+XMU*V*COSU+G3*Gt)*(  .25*SINUMQ) 

HVV={.01)*2.*B1*FI1*G5*(1.-«25*COSUMQ) 

HVU  =  -XLA.M*SINU+XMU*COSU+G5*.01*FI1*2.*B1*V*( .2  5*SINUMQ) 

HUU=-XLAM*V*C0SU-XMU*V*SINU+G3*G5*( .2  5*COSUMO) 

HW=G3*( l.-.2  5*(COSUMQ)  ) *G4+F I  2* ( -2 .* ( W/Wl ) * ( 1./W1) )/(G7*G7) 

HWU=G3*( .25*(SINUMQ) )*G4 

HWV=G4*(  .01  )*FI1*(2.*B1*V)*(1.-.2  5*(C0SUMQ)  ) 

HWW=   G3*( ( l.-.2  5*COSUMQ)*(Pl+P2+P3+P4) )+P9 

P5-HVV*HWW-HWV*HWV 

P6=(-HV)*HWW+HW*HWV 

P7-HVV* (-HW ) +HWV*HV 

DETER=HUU*P5-HVU* ( HVU*HWW-HWU*HWV ) +HWU* ( HVU*HWV-HWU*HVV ) 

FUM1=(-HU)*P5-HVU*P6+HWU*( (-HV ) *HWV+HW*HVV ) 

FUM2=HUU*P6+HU*  (  HVU*HWW-HWU*HWV  )  +HWU*  (  HVU*  (  -HW  )  +HWU*HV  ) 

FUM3=HUU*P7-HVU* ( HVU* ( -HW ) +HWU*HV ) -HU* ( HVU*HWV-HWU*HW ) 

CONTINUE 

RETURN 

END 
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SUBROUTINE    BOUNDW ( U» V » W»XLAM»XMU»FI 1 »A1 ,B1 »C1 »D1 , WO » FI 2 »W1 »Q »WMAX » 
+HW ) 
W=WMAX  ] 

COSU=COSF(U)  1 

SIi\U  =  SINF(U)  | 

COSUMQ=COSF(U-Q) 
SINUMQ=SINF(U-Q) 
Gl=Cl+Dl*(W-WO)*(W-WO) 
G2=l.+Dl*(W-WO)*(W-WO) 
G3=(.01)*FI1*(A1+B1*V*V) 

G4=(-G1*(2.*D1*(W-W0) )+G2*(2.*Dl*(W-WO) ) )/(G2*G2> 
G5=(Ci+Dl*(W-WO)*(W-WO) )/(l.+Dl*(W-WO)*(W-WO) ) 

G7=1.+(W/W1)**2  I 

HV=XLAM*COSU+XMU*SINU+( .01 ) *F 1 1* ( 2 .*B1*V ) *G5*( 1 .-.2  5*COSUMQ) 
HU=^XLAM*V*SINU+XMU*V*COSU+G3*G5*( .2  5*SINUMQ) 
HVV=( •01)*2.*B1*FI l*G5*(l.-.2  5*COSUMQ) 

HVU=-XLAM*SINU+XMU*COSU+G5*.01*FI1*2.*B1*V*( .25*SINUMQ) 

HUU=-XLAM*V*C0SU-XMU*V*SINU+G3*G5*(  .25*COSUMQ)  Bj 

HW=G3*( i.-.2  5*(COSUMQ) ) *G4+FI 2* ( -2 .* ( W/Wl ) *( 1./W1) )/(G7*G7) 
DEN=HVU*HVU-HUU*HVV 

GUM1=-HV*HVU+HU*HVV  | 

GUM2=HVU*(-HU)+HV*HUU  I 

DO  981  1=1,7  | 

DU=GUM1/DEN 
DV=GUM2/DEN 
U=U+DU 
V=V+DV 
COSU=COSF(U) 
SINU=SINF(U) 
COSUMQ=COSF(U-Q) 
SINUMQ=SINF(U-Q) 

G1=C1+D1*(W-W0)*(W-W0)  | 

G2=1.+D1*{W-WO)*(W-W01 

G3=( •01)*FI1*( A1+B1*V*V)  | 

G4=(-G1*(2.*D1*(W-W0)  )+G2*(2.*Dl*(W-WO)  )  }/(G2*G2) 
G5=(C1+D1*(W-W0)*(W-W0) )/( l.+Dl*(W-WO)*(W-WO) ) 


• 
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HV=XLAM*COSU+XMU*SINU+( .01 )*Fl 1* ( 2 .*B1*V ) *G5*( 1 .-.25*COSUMQ ) 
HU=-XLAM#V*SINU+XMU*V*C0SU+G3*G5*( .25*SINUMQ ) 
HVV=( .01)*2«*B1*FI1*G5*< 1 .-. 2  5*COSUMQ ) 

HVU=-XLAM*SlNU+XMU*COSU+G5*.0l*FIl*2.*Bl*V*( •2  5*SINUMQ) 
HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .2  5*COSUMQ) 

:n-hvu*hvu-huu*hvv 

GUM1=-HV*HVU+HU*HVV 

GUM2=HVU*  (  -HU  )  +HV*HUU 
981  CONTINUE 
RETURN 
END 


SU3R0UTINE  SOUNDV ( U, V > W,XLAM»XMU» FI 1 ♦ Al ,B1 »C1 »D1 . WO. FI 2 »Q» VMAX. 
+  W1 ) 
V=VMAX 

COSU=COSF(U) 
SINU=SINF(U) 
COSUMQ=COSF(U-Q) 
SINUMQ=SINF(U-Q) 
G1=C1+D1*(W-W0)*(W-W0) 
G2=l.+Dl*(W-WO)*(W-WO) 
G3=( .01)*FI1*(A1+B1*V*V) 

G4={-G1*(2.*D1*(W-W0) )+G2*(2.*Dl*(W-WO) ) )/(G2*G2) 
G5=(C1+D1*(W-W0)*(W-W0) )/( l.+Dl*(W-WO)*(W-WO) ) 
G7=1.+W*W/(W1*W1) 

Pl=( (l.+Dl*(W-WO)*(W-WO) )*2.*D.1+4.*D1*(W-W0)*D1*(W-W0) )/(G2*G2) 
P2=(-(  8.*Dl*(V/-WO)*(W-WO)  )*Di)/(G2*G2) 

P3=     (-{Cl+Dl*(W-WO)*(W-WO> )*2.#D1-4.*D1*(W-W0)*D1*(W-W0> )/(G2*G; 
P4=     (G1*8.*D1*(W-W0)*D1*(W-W0) J / ( G2*G2*G2 ) 
P9=FI2*(-2«/(Wl*Wl*G7*G7)+8.*W*W/(Wl*Wl*Wl*Wl*G7*G7*G7) ) 
HU=-XLAM*V*SINU+XMU*V*COSU+G3*G5*( •25*SINUMQ) 
HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .2  5*COSUMQ) 
HW=G3*(  l.-.2  5*(COSUMQ)  )*G4+FI2* (-2 •*( W/Wl )*<  1./W1  )  )/(G7*G7) 
HWU=G3*(.25*(SINUMQ)  )*G4 
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HWW=   G3*( ( l.-.25*COSUMQ) *(P1+P2+P3+P4) )+P9 

DEN=HWU*HWU-HUU*HWW 

GUM1=-HW*HWU+HU*HWW 

GUM2=-HU*HWU+HW*HUU 

DO  981  1=1,7 

DU=GUM1/DEN 

DW=GUM2/DEN 

U=U+DU 

W=W+DW 

cosu=cosF(in 

SINU»SINF(U) 

COSUMQ=COSF(U-Q) 
SINUMQ=SINF(U-Q) 
Gl=Cl+Dl*(W-WO)*(W-WO) 
G2=l.+Dl*(W-WO)*(W-WO) 
G3=( .01)*FI1*(A1+B1*V*V) 

G4=(-G1*(2.*D1*(W-W0) )+G2*(2.*Dl*(W-WO) ) )/(G2*G2) 
G5=<C1+D1*(W-W0)*(W-W0)  }/  (  l.+Dl*(W-WO)*(W-WO)  ) 
G7=1.+W*W/(W1*W1) 

Pl=( (l.+Dl*(W-WO)*(W-WO) )*2.*D1+4.*D1*(W-W0)*D1*(W-W0) )/(G2*G2) 
P2=(-(  8«*Dl*(W-WO)*.(W-WO)  )*Di)  /(G2*G2) 

P3  =  (-(Cl+Dl*(W-WO)*(W-WO) )*2.*Dl-4.*Dl*(W-WO)*Dl*(W-WO) )/(G2*G2) 
P4=  (Gl*8«*Dl*(W-WO)*Dl*(W-WO) ) /(G2*G2*G2) 
P9=FI2*(-2./(Wl*Wl*G7*G7)+8.*W*W/(Wl*Wl*Wl*Wl*G7*G7*G7> ) 
HU=-XLAM*V*S I NU+XMU*V*COSU+G3*G5* ( .  25*S I NUMQ ) 
HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .25*COSUMQ) 
HW=G3*( l.-.2  5*(C0SUMQ) ) *G4+FI 2* < -2 .* ( W/Wl ) *( 1./W1) )/(G7*G7) 
HWU=G3*< .2  5*(SINUMQ) )  *G4 

HWW=   G3*( ( l.-.25*COSUMQ)*(Pl+P2+P3+P4) )+P9 
DEN=HWU*HWU-HUU*HWW 
GUM1=-HW*HWU+HU*HWW 
GUM2=-HU*HWU+HW*HUU 
981  CONTINUE 
RETURN 
END 
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SUBROUTINE  SEARCH ( Al »B1 »C1 »D1 , FI 1 » FI 2 »XLAM .XMU.COSU.U.V ,W.S INU» 
+WO  ,  G 1 »  G  2  » W 1 » COSUMQ  ) 
DIMENSION  VINT(10j,WINT(10),H(50,50) 
F1V=(A1+81*V*V)*< .01)*(1.-.2  5*COSUMO)*G1/G2 
F2V»1./(1#+(W*W/(W1*W1 ) ) ) 

HAM     =XLAM*V*C0SU+XMU*V*.SINU  +  F1V*FI1  +  F2V*FI2 
VINT(0)=0.0 
WINT(0)=0.0 
DO  131  M=1.5 
VINT(M)=VINT(M-l)+3.0 
DO  122  N=l,5 
WINT(N)=WlNT(N-l)+200.0 
G 1  =  C 1  +  D 1* ( W I N T ( N ) -WO ) * ( W I  NT ( N ) -WO  ) 
G2  =  1.+D1*(WINT(N) -WO ) * ( W I  NT ( N ) -WO ) 

F1V=(A1+S1*VINT(M-)*VINT(M)  )*(  .01  >  *  (  1  •-•25*COSUMQ)  *G1/G2 
F2 V  =  l. /( l.+W I  NT (N)*WI NT (N)/(W1*W1)  ) 
H(M»N)=XLAM*VINT(M)*C0SU+XMU*VINT(M)*SINU+F1V*FI1+F2V*FI2 

122  CONTINUE 
131  CONTINUE 

HAMSRCH=H( 1*1) 
VSRCH=VINT( 1) 
WSRCH=WINT( 1) 
DO  121  M=l»5 
DO  123  N=1.5 

IF  (HAMSRCH-H(M.N) )  123.123*111 
111  HAMSRCH=H(M»N) 
VSRCH=VINT(M) 
WSRCH=WINT(N) 

123  CONTINUE 
121  CONTINUE 

IF  (HAMSRCH-HAM)  1.2*2' 

1  V=VSRCH 
W=WSRCH 

2  RETURN 
END 
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SUBROUTINE  WORSI  ( U» V, W»XLAM»XMU»F I 1 » Al ,B1 »C1 tDl »WO» FI2 tWl »0) 

W=500. 

SINU=SINF(U5 

COSUMQ=COSF(U-Q) 

SINUMQ=SINF(U-Q) 

Gl =C i+D 1* ( W-WO ) * ( W-WO ) 

G2 = 1. +D 1* ( W-WO )*( W-WO) 

G3=U0l)*FIl*(Al+Bl*V*V) 

G4=(-G1*(2«*D1*{W-W0) )+G2*(2.*Dl*(W-WO) ) )/(G2*G2) 

G5=(C1+D1*( W-WO)* (W-WO) )/( I. +01* (W-WO)* (W-WO) ) 

G7  =  l.+W/..0 

Pl  =  J (1.+D1* (W-WO)* (W-WO) ) *2. *0 1+4 • *D 1* ( W-WO )*D1*( W-WO) )/(G2*G2) 

P2=(-( 8.*Di* (W-WO) *( W-WO) )*D1) /(G2*G2) 

P3=  (-( CI +D1*( W-WO)* (W-WO) ) *2. *D1-4.*D1*( W-WO )*D1*( W-WO) )/(G2*G2) 

P4=  (G1*8.*D1*(W-W0)*D1*(W-W0) ) /(G2*G2*G2) 

P9=FI2*( ( (-G7*( .1/W0)+G6*(1./W0) )*2.*(1./W0) )/(G7*G7*G7) ) 

HV  =  XLAM*COSU+XMU*SINU+(.01 ) *F 1 1* { 2 .*B1*V ) *G5* ( 1 .-.25*COSUMQ ) 

HU=-XLAM*V*SINU+XMU*V*COSU+G3*G5*< .2  5*SINUMQ) 

HVV=(  .01)*2.*B1*FI  1*G5*(  1  .-.25*C0SUMQ  ) 

HVU=-XLAM*SINU+XMU*COSU+G5*.0l*FIl*2.*Bl*V*( .25*SINUMQ) 

HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .25*C0SUMQ) 

HW=G3*< l.-.25*(C0SUMQ) )*G4+(Fl?*( (G7*( .1/WO) )-G6*(l./W0) ) )/(G7*G7) 

HWW=   G3*(  (  l.-.25*C0SUMQ)*>(Pl  +  P2+P3+P4)  )+P9 

DEN=HVU*HVU-HUU*HVV 

GUM1=-HV*HVU+HU*HVV 

GUM2=HVU* ( -HU ) +HV*HUU 

DO  981  1=1,7 

DU=GUM1'/DEN 

DV=GUV;2/DEN 

U=U+DU 

V=V+DV 

COSU=COSF(U) 

SINU=SINF(U) 
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COSUMQ=COSF(U-G) 
SINUMQ=SINF(U-Q) 
G 1 =C 1+0 1* ( W-WO ) * ( W-WO ) 
G2=l »+D 1*< W-WO )*( W-WO) 
G3=(.01)*FI 1*(A1+B1*V*V) 

G4=(-G1*(2.*D1*(W-W0) )+G2*(2.*Dl*(W-WO) ) )/(G2*G2) 
G5=(C1+D1*( W-WO)* (W-WO) )/(l.+Dl*(W-WO)*(W-WO) ) 
F1V«(A1+B1*V*V)*< .01 )*(1.-.2  5*C0SUMQ)*G1/G2 
F2V=(l.+.l*W/WO)/( l.+W/WO) 

HV=XLAM*COSU+XMU*SINU+( .01 )*FI 1*(2 .*31*V)*G5*( l.-«25*COSUMQ) 
HU=-XLAM*V*SINU+XMU*V*COSU+G3*G5*( .25*SINUMQ) 
HVV=( .01)*2.*B1*FI1*G5*( 1 .-. 2  5*COSUMQ ) 

HVU=-XLAM*SlNU+XMU*COSU+G5*.0l*FIl*2.*Bl*V*( .2  5*SINUMQ) 
HUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .2  5*COSUMQ) 
DEN=HVU*HVU-HUU*HVV 
GUM1=-HV*HVU+HU*HVV 
GUM2*HVU* ( -HU ) +HV*HUU 
981  CONTINUE 
RETURN 
END 
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Appendix  II 

In  this  appendix  additional  information  about  the 
paths  described  in  section  7  and  the  path  discussed  in 
section  8  \nll   be  given.  For  the  first  three  paths 
the  position  (x,y)  and  the  controls  u,v,w  for  the  sub- 
marine will  be  given  at  twelve -hour  intervals.  This 
information  will  be  given  for  Path  IV  and  in  addition 
the  submarine's  location  and  controls  will  be  included 
for  the  time  steps  immediately  proceeding  and  following 
the  corner. 

Path  I 

Time     x       y_       u       v       w 

0.0  0.0  0.0  -  0.7 

12.0  104.5  -  74.6  -  0.6 

24.0  214.8  -■140,4  -  0.5 

36.0  330.0  -  197.2  -  0.4 

48.0  449.2  -  244.8  -  0.3 

60.0  57L5  -  283.1  -  0.3 

72.0  696.1  -  312.5  -  0.2 

84.0  822.3  -  333.0  ■  -  0,1 

96.0  949.3  -  345.1  -  0.1 


10.7 

207.8 

10.7 

207.8 

10.7 

207.8 

10.7 

207.8 

10.7 

207.7 

10.7 

207.5 

10.7 

207.3 

10.6 

207.1 

10.6 

206.8 
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U  V  w 

0.0  10.6  206.5 

0.1  10.6  206.2 

0.1  10.5  205.9 

0.2  10.5  205.5 

0.2  10.5  205.2 

0.3  10.5  204.9 

0.3  10.4  204.6 

0.3  10.4  .  204.3 

0.4  10.4  204.0 

0.4  10.4  203.7 

0.4  10.3  203.5 

0.5  11.1  203.7 

Path  II 

Time      x       £       u  v  w 

0.0     0.0     0.0   -  0.6  10.1  1000.0 

12.0    102.4   -  66.1   -  0.5  10.2  1000.0 

24.0    210.1  -  124.3   -  0.5  10.2  1000.0 

36.0    322.4  -  174.3   -  0.4  10.3  1000.0 

46.0    438.6  -  216.3  '  -  0.3  10.3  1000.0 

60.0    558.0  -  250.4   -  0.2  10.4  1000.0 

72.0    680.0  -  276.7   -  0.2  10.4  1000.0 


108.0 

1076.5 

-  349.2 

120.0 

1203.5 

-  345.7 

132.0 

1329.8 

-  335.1 

144.0 

1455.0 

-  318.0 

156.0 

1578.9 

-  294.8 

168.0 

1701.3 

-  266.0 

180.0 

1822.0 

-  232.I 

192.0 

1941.0 

-  193.6 

204.0 

2058.2 

-  150.8 

216.0 

2173.6 

-  104.2 

228.0 

2287.2 

-  5^.1 

240.0 

2400.0 

0.0 
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Tiiaa 


x 


u 


w 


"" 

" 

84.0 

804.0 

-  295.6 

-  0.1 

10.5 

1000.0 

96.0 

929.4 

-  307.4 

-  0.1 

10.5 

1000.0 

1C6.0 

1055.7 

-  312.3 

0.0 

10.5 

1000.0 

120.0 

1182.4 

-  310.9 

0.0 

10.6 

1000.0 

132.0 

1309.2 

-  303.7 

0.1 

10.6 

1000.0 

144.0 

1435.7 

-  290.9 

0.1 

10.6 

1000.0 

156.0 

1561.7 

-  273.2 

0.2 

10.6 

1000.0 

168.0 

1686.9 

-  250.8 

0.2 

10.6 

1000.0 

180.0 

1811.2 

-  224.1 

0.2 

10.6 

1000.0 

192.0 

1934.5 

-  193.7 

0.3 

10.6 

1000.0 

204.0 

2056.6 

-159.3 

0.3 

10.6 

1000.0 

216.0 

21?7.6 

-  122.7 

0.3 

10.5 

1000.0 

228.0 

2290.2 

-  64.0 

0.5 

10.6 

705.2 

240.0 

2400.0 

0.0 

Path 

0.5 

III 

10.7 

650.1 

Tima 


u 


w 


0.0 

0.0 

0.0 

-'0.7 

10.8 

527.4 

12.0 

106.1 

-  75.4 

-  0.6 

10.9 

526.4 

24.0 

218.1 

-  141.9 

-  0.5 

10.9 

526.6 

36.0 

334.9 

-  199.1 

-  0.4 

10.8 

528.1 

48.0 

455.8 

-  246.9 

-  0.3 

10.8 

530.7 
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Tinio  x  y  u 


60.0 

579.6 

-  285.4 

-  0.3 

10.8 

5>.6 

72.0 

705.5 

-  314.6 

-  0.2 

10.8 

539.6 

84.0 

832.8 

-  33^.8 

-  0.1 

10.7 

5^5.9 

96.0 

960.6 

-  346.5 

-  0.1 

10.7 

553.3 

108.0 

IO88.3 

-  350.1 

0.0 

10.6 

561.9 

120.0 

1215.4 

-  346.1 

0.1 

10.6 

571.7 

132.0 

1341.5 

-  335.0 

0.1 

10.5 

582.7 

144.0 

1466.3 

-  317.4 

0.2 

10.5 

594.9 

136.0 

1589.5 

-  293.9 

0.2 

10.4 

608. 3 

168.0 

1710.9 

-  264.8 

0.3 

10.4 

622.9 

180.0 

I830.5 

-  230.7 

0.3 

10.3 

638.7 

192.0 

1948.1 

-  192.1 

0.3 

10.3 

655.8 

204.0 

2063.8 

-  149.4 

0.4 

10.3 

674.2 
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-  103.0 
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10.2 

693.9 
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2289.6 

-  53.2 
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10.2 

715.0 

240.0 

2400.1 

0.0 

0.5 

10.4 

677.9 

Path  IV 

Tine  x  £  u  v    *  w 

0.0  0.0  0.0       -  1.1  11.0        200.8 

12.0  70..2     -  111.9      -  1.0  11.0        200.8 
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24.0 
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-  218.0 

-  0.9 

11.1 

200.8 

36.0 

237.7 

-  317.5 

-  0.8 

11.1 

200.8 

48.0 

334.1 

-  409.4 

-  0.7 

11.1 

200.8 

60.0 

437.9 

-  493.1 

-  0.6 

11.1 

200.8 

72.0 

543.4 

-  568.0 

-  0.6 

11.1 

200.9 

84.0 

664.6 

-  633.6 

-  0.5 

11.1 

200.8 

96.0 

785.7 

-  689.6 

-  0.4 

11.1 

200.8 

108.0 

910.5 

-  736.0 

-  0.3 

11.1 

200.8 

120.0 

1038.3 

-  772.7 

-  0.2 

11.1 

200.8 

132.0 

1168.3 

-  799.8 

-  0.2 

11.1 

200.8 

144.0 

1302.0 

-  816.3 

-  0.1 

11.5 

201.1 

156.0 

1453.5 

-  811.1 

0.2 

15.0 

202.9 

158.0 

1482.9 

-  804.7 

0.3 

15.0 

204.5 

160.0. 

15H.5 

-  796.0 

0.3 

15.0 

207.0 

162.0 

1539.4 

-  784.8 

0.4 

15.0 

211.4 

164.0 

1565.1 

-  771.6 

0.5 

15.0 

1000.0 

166.0 

159L5 

-  757.2 

0.5 

15.0 

1000.0 

168.0 

1617.3 

-  742.0 

0.6 

15.0 

1000.0 

180.0 

1759.6 

-  632.3 

0.7 

15.0 

1000.0 

192.0 

1889.7 

-  507.8 

0.8 

15.0 

1000.0 

204.0 

2017.5 

-  381.1 

0.8 

15.0 

1000.0 

216.0 

2145.1 

-  254.1 

0.8 

15.0 

1000.0 

228.0 

2272.6 

-  127.0 

0.8 

15.0 

1000.0 

240.0 

2400.0 

0.1 

0.8 

15.0 

1000.0 
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